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Abstract. We present a rigorous derivation of a real space Full- Potential Multiple- 
+^ Scattering- Theory (FP-MST) that is free from the drawbacks that up to now have 

S impaired its development (in particular the need to expand cell shape functions in 
spherical harmonics and rectangular matrices), valid both for continuum and bound 
states, under conditions for space-partitioning that are not excessively restrictive and 
easily implemented. In this connection we give a new scheme to generate local basis 
functions for the truncated potential cells that is simple, fast, efficient, valid for any 
shape of the cell and reduces to the minimum the number of spherical harmonics in 
the expansion of the scattering wave function. The method also avoids the need for 
saturating 'internal sums' due to the re-expansion of the spherical Hankel functions 
around another point in space (usually another cell center). Thus this approach, 
' ^S^ provides a straightforward extension of MST in the Muffin-Tin (MT) approximation, 

with only one truncation parameter given by the classical relation Z max = kRf,, where 
k is the electron wave vector (either in the excited or ground state of the system 
under consideration) and Rb the radius of the bounding sphere of the scattering cell. 
Moreover, the scattering path operator of the theory can be found in terms of an 
absolutely convergent procedure in the Z max — > oo limit. Consequently, this feature 
provides a firm ground to the use of FP-MST as a viable method for electronic structure 
calculations and makes possible the computation of x-ray spectroscopies, notably 
$_i photo-electron diffraction, absorption and anomalous scattering among others, with 

the ease and versatility of the corresponding MT theory. Some numerical applications 
of the theory are presented, both for continuum and bound states. 



o 
o 



PACS numbers: 78.70Dm, 61.05.jd 
1. Introduction 

At its most basic, Multiple scattering Theory (MST) is a technique for solving a linear 
partial differential equation over a region of space with certain boundary conditions. It 
is implemented by dividing the space into non-overlapping domains (cells), solving the 
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differential equation separately in each of the cells and then assembling together the 
partial solutions into a global solution that is continuous and smooth across the whole 
region and satisfies the given boundary conditions. 

As such MST has been applied to the solution of many problems drawn both 
from classic as well as quantum physics, ranging from the study of membranes and 
electromagnetism to the quantum-mechanical wave equation. In quantum mechanics 
it has been widely used to solve the Schrodinger equation (SE) ( or the associated 
Lippmann-Schwinger equation (LSE)) both for scattering and bound states. It was 
proposed originally by Korringa and by Kohn and Rostoker (KKR) as a convenient 
method for calculating the electronic structure of solids [TJ [2] and was later extended 
to polyatomic molecules by Slater and Johnson [3]. A characteristic feature of the 
method is the complete separation between the potential aspect of the material under 
study, embodied in the cell scattering power, from the structural aspect of the problem, 
reflecting the geometrical position of the atoms in space. 

Applications of the KKR method were first made within the so-called muffin- 
tin (MT) approximation for the potential. In this approximation the potential is 
confined within non-overlapping spheres, where it is spherically symmetrized, and takes 
a constant value in the interstitial region. Moreover, although spherical symmetry is not 
formally necessary, the condition that the bounding spheres do not overlap was thought 
to be necessary for the validity of the theory. Despite this approximation the method 
is complicated and demanding from a numerical point of view and as a band-structure 
method it was therefore superseded by more efficient linearized methods, such as the 
linearized muffin-tin-orbital method (LMTO) [I] and the linearized augmented-plane- 
wave method (LAPW). 

Full-potential versions of these band methods have also been introduced in recent 
years. However, none of these methods can match the power and versatility of a full- 
potential method based on the formalism of MST, either in terms of providing a complete 
solution of the SE or in the range of problems that could be treated. In particular, none 
of these methods leads easily to the construction of the Green's Function (GF) which is 
indispensable in the study of a number of properties of many physical systems. 

Due to these reasons, in the last two decades, the KKR method has experienced 
a revival in the framework of the Green's function method (KKR-GF). Indeed, due 
to the introduction of the complex energy integration, it was found that the method 
is well suited for ground-state calculations, with an efficiency comparable to typical 
diagonalization methods. An host of problems became in this way tractable, ranging 
from solids with reduced symmetry (like e.g. isolated impurities in ordered crystal, 
surfaces, interfaces, layered systems, etc..) to randomly disordered alloys in the coherent 
potential approximation (CPA). 

At the same time it soon became clear that the MT approximation was not adequate 
to the treatment of systems with reduced symmetry or for the calculation of lattice forces 
and relaxation. In order to deal with these problems a number of groups developed a 
full potential (FP) KKR-GF method, obtaining very good results, comparable with 
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full-potential LAPW method (FLAPW), for what concerns total energy calculations, 
lattice forces, relaxation around an impurity, ( [61 13 [10] and Refs. therein). 
Due to their method of generating the single site solutions and the cell t-matrix, the 
additional numerical effort required for the implementation of the FP-MS scheme scales 
only linearly with the number of non-equivalent atoms and is not significantly greater 
than in the MT case. 

In this development the authors took an empirical attitude toward some 
fundamental problems related to the extension of MST to the full-potential case, like 
the strongly debated question of the /-convergence of the theory or the need to converge 
"internal" sums arising from the re-expansion of the free Green's function around two 
sites, which entails the unwanted feature of the introduction of rectangular matrices 
into the theory. [TT] Without getting involved into ab initio questions, they just use 
square matrices for the structural Green's Function G 7 l r l,(E) needed to calculate the 
Green's Function of the system (see e.g. Eqs. (6) and (9) in Ref. [9]) and truncate the 
/-expansion to Z max = 3 or 4, obtaining in this way the same accuracy as the FLAPW 
method. 

Some observations are in order at this point. First, the FP method in the framework 
of MST has been initially developed only for periodic systems in two or three dimensions 
and for states below the Fermi level. To our knowledge, its extension to treat bound and 
continuum states of polyatomic molecules and in general real space applications of the 
method have progressed very slowly and have been scarce. Secondly, the generation of 
the local solutions of the SE with truncated cells in the FP extension of the MST has up 
to now involved the expansion of the cell shape function in spherical harmonics, which 
might create convergence problems, as discussed below. Thirdly, the FP extension of 
MST has generated a lot of controversies that have gone on for more than thirty years 
|12j . Some of the problems have found a solution and we refer the reader to the book of 
Gonis and Butler [13] for a comprehensive review of the state of the art in this field (in 
particular see their chapter 6). However, questions like the /-convergence of the theory 
or the use of square matrices are still matter of debate and some rigorous answer should 
be given to them. 

As mentioned above, applications to states well above the Fermi energy, as required 
in the simulations of x-ray spectroscopies, like absorption, photo-emission, anomalous 
scattering, etc., have been scarce. In the words of Ref. [13], "the feeling that one 
should calculate the "near-field-corrections" (NFC), coupled with the need to solve a 
fairly complicated system of coupled differential equation to determine the local (cell) 
solutions (based on the phase function method) has contributed greatly to the slow 
development of a FP method based on MST". It was only after it was realized that 
NFC are not necessary and a new method to generate local solutions was found that 
progress became faster, at least in the calculation of the electronic structure of solids. 
( 0, El O HO] and Refs. therein) The only remaining drawback was and is the 
truncation of the potential at the cell boundary which is still performed via a shape 
function expanded in spherical harmonics. Added to this there is the feeling that one 
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should still converge the "internal" sums leading to the use of rectangular matrices in the 
angular momentum (AM) indexes, although this last step is sometimes ignored without 
justifications. Last, but not least, the question of the /-convergence of the theory remains 
unsettled. 

For all these reasons FP codes based on MST for the calculation of x-ray 
spectroscopies are not very numerous. We mention here the work by Huhne and Ebert 
[13] on the calculation of x-ray absorption spectra using the FP spin-polarized relativistic 
MST and that of Ankudinov and Rehr [15] in the scalar relativistic approximation. 
These authors use the potential shape function to generate the local basis functions 
which are at the heart of MST. The expansion of the shape function and the cell potential 
in spherical harmonics leads to a high number of spherical components in the coupled 
radial equations that becomes progressively cumbersome to handle and time consuming 
with increasing energy and in absence of symmetry. This feature might also be at the 
origin of another problem related with the saturation of "internal" sums in the MSE 
[13] . as discussed later in this paper. Moreover no critical discussion is devoted in their 
work to the /-convergence problems of MST or the use of square matrices in the theory. 

Another code based on a version of the MST that uses non overlapping spherical 
cells and treats the interstitial potential in the Born approximation is that of Foulis et 
al. [TBI [T7] This method however treats in an approximate way the potential in the 
interstitial region and moreover looses one of the major advantages of the MST, namely 
the separation between dynamics and geometry in the solution of the scattering problem. 
Foulis [18] is now developing an exact FP-MS scheme based on distorted waves in the 
interstitial region that seem to be promising, but its numerical implementation is still 
to come. 

There are other codes that simulates x-ray spectroscopies and are not based on 
MST: that of Joly [TH] is based on the discretization of the Laplacian in three dimensions 
(finite-difference method (FDM)), where the SE is solved in a discretized form on a three- 
dimensional grid, the values of the scattering wave-function being the unknowns. This 
method is however limited to cluster sizes of the order of 20 atoms (without symmetry), 
due to the high memory requirement when the number of mesh points increases with 
the dimensions of the cluster. Finally a method based on the pseudo-potential theory 
to calculate x-ray absorption is worth mentioning. [20] It can easily cope with clusters 
of many atoms (300 and more) with a computational effort that scales linearly with the 
number of atoms. One of its drawback is its little physical transparency and the fact 
that it has been applied only to calculate x-ray absorption spectra. Also, relaxation 
around the core hole must be taken into account by super-cell calculations and there is 
little flexibility to deal with energy-dependent complex potentials. 

The purpose of the present paper is the rigorous derivation of a real space FP-MST, 
valid both for continuum and bound states, that is free from the drawbacks hinted to 
above, in particular the need to use cell shape functions and rectangular matrices, 
under conditions for space-partitioning that are not excessively restrictive and easily 
implemented (see beginning of Section [3]). In connection with this we shall present a 
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new scheme to generate local basis functions for the truncated potential cells that is 
simple, fast, efficient, valid for any shape of the cell and reduces to the minimum the 
number of spherical harmonics in the expansion of the scattering wave function. Finally 
we shall also address the problem of the /-convergence of the theory, giving a positive 
answer to this debated question. 

Even though this work is primarily motivated by applications in spectroscopy, it 
will be clear from the context that bound states can be treated as well. Actually the 
method can also work for complex energy values, so that one can take advantage of the 
fact that the solution of the Schrodinger equation is analytical in the energy plane, as 
is the associated Green's function, except for cuts and poles on the real axis. Therefore 
spectroscopy is only one regime of applications. 

Section [2] of this paper presents the new scheme to generate local basis functions and 
tests it against known solutions for potentials cells with and without shape truncation. 
Section [3] provides a new derivation of the FP-MST that allows us to work with square 
matrices for the phase functions Sll' and Ell> and for the cell Tll> matrix with only 
one truncation parameter, contrary to the present accepted view. [13] Due to their 
importance in the theory, various equivalent forms for the Green's function are presented 
in this scheme. This latter is extended to the calculation of bound states of polyatomic 
molecules and tested against the known eigenvalues of the hydrogen molecular ion. 
Section [4] discusses the strongly debated problem of the /-convergence of the theory and 
provides a truncation procedure that converges absolutely in the Z max — > oo limit. 

Section [5] reports one additional application of the present FP-MS theory besides 
those already presented in Ref.s [21] and |22j . namely the calculation of the absorption 
cross section in the case of linear molecules {Br 2 diatomic molecule), where the 
improvement of over the MT approximation is quite dramatic. Moreover, with an eye 
to using the theory to study the performance of model optical potentials, Section [5] also 
presents a preliminary application of the non-MT (NMT) approach to the study of the 
relative performance of the Hedin-Lundqvist (HL) and the Dirac-Hara (DH) potentials 
in the case of a transition metal. Finally Section[6]presents the conclusions of the present 
work. A preliminary and partial account of this latter has been presented in Ref.s [21] 
and [22] • 

2. Local Basis functions for single truncated potential cells 

A characteristic feature of MST is that it does not rely on a finite basis set for the 
expansion of the global wave function inside each cell as all other methods of electronic 
structure calculations do. Instead it relies on expanding the global solution in terms 
of local solutions of the Schrodinger equation at the energy of interest, which can be 
regarded as an optimally small, energy adapted basis set [9]. Therefore it is essential 
for the practical implementation of the theory to devise an efficient numerical method 
to generate them. We shall consider Williams and Morgan (WM) basis functions 
[23J which inside each cell are local solutions of the SE and behave at the origin as Jl{ v ) 
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for r — > 0. Throughout the paper we shall use real spherical harmonics and shall put for 
short J L (v-k) = ji(kr)Y L (r), N L (r-k) = ni (kr)Y L (r) and H^(t; k) = -ikhf (kr)Y L (v) , 
where ji,n>i,hi denote respectively spherical Bessel, Neumann and Hankel functions of 
order I. The truncated cell potential V(r, r) is defined to coincide with the global 
system potential inside the cell and to be equal to zero (or to a constant) outside. As 
mentioned in the introduction we want to avoid the expansion of the truncated cell 
shape function (or equivalently of the truncated potential) in spherical harmonics due 
to convergence problems. However we observe that, even if the potential has a step, 
the wave function and its first derivative are continuous, so that its angular momentum 
expansion is well behaved and even converges uniformly in r. [24J Therefore we can 
safely write $i(r) = YIl' Rl'l{t)Yl'{y) and this expression can be integrated term by 
term under integral sign. 



2.1. Three-dimensional Numerov method 



In order to generate the basis functions we write the SE in polar coordinates for the 
function Pl{ v ) — 



P L (r,r) = -~L 2 P L (r,v) 



(1) 



where L 2 is the angular momentum operator, whose action on Pl(t, r) can be calculated 
as: 



L 2 P L (r,i)=J2ni' + l)rRu L {r)Y LI {v) 



(2) 



Equation (JI]) in the variable r looks like a second order equation with an 
inhomogeneous term. Accordingly we use Numerov's method to solve it. As is well 
known, putting /y = Pt(rj, ij) and dropping for simplicity the index L, the associated 
three point recursion relation is 

h 6 



A.+i,jfi+i,j - Bijfij + Ai-ijfi-ij — gij - ^^fi,j 



(3) 



where, 



.4 



B 



1,3 



V 



'■J 



9i,j 

%3 



1 ^ 

2 + —v tJ = 12 - 10A itj 

v{r i} Tj) - e 

h 2 

— fe+ij + lOgy + qi-ij] 
\Y, l '( l ' + 1 ^ R L'L(n)Y L '(r :i 



(4) 



Here % is an index of radial mesh and j an index of angular points on a Lebedev 
surface grid. [25] Obviously rji?L^(rj) = "%2jWjPi,(ri,Tj)YLr(Tj), where Wj is the 
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weight function for angular integration associated with the chosen grid. The number 
of surface points N Leb is given by N Leb « (2Z max + l) 2 /3 as a function of the maximum 
angular momentum used [26] , taking into account that one integrates the product of two 
spherical harmonics. As it is, we cannot use Eq. ^ to find fi+ij by iteration, from the 
knowledge of and fi-ij at all the angular points, since the " inhomogeneous" term 
qi+ij is not expressible in terms of fi+ij due to the last line of Eq. Q and is calculated 
at the radial mesh point i + 1. 

We first eliminate this point from the expression of gij, observing that 

h 2 

9i,j ~ — 



12 
12 



[q i+ i,j + lOqi + qi-i,j] 

Qi+l,j ~ ^Qi + Qi-l,j ,2 

h 2 



12q 



The second order central difference is given by 



q i+1 - 2qi + q 



so that 



12 



i-l 



1 I :J 



h 2 q'l + ^qT + i^q? + 



h 4 



12 



360 



20160 



qr + 



h 2 

12' 



(5) 



(6) 



(?) 



omitting errors of order h % and higher. 

Now for the second derivative q"^ we use the backward formula ^27\ 



qu - 2qi-i,j + q, 



+ k; 



7h 2 



h 2 ' 12 

to avoid the contribution of the point % + 1. Inserting Eq. ^ into Eq. ^ 



h 2 



h 5 



h 6 



9 



(8) 



(9) 



which is the formula we wanted to arrive at. Therefore our modified Numerov procedure 
becomes: 

h 5 



Ai+ijfi+ij — Bijfij + Aj_i/j_x,j — g%j + — q 



(10) 



where, 



Ay — 1 — 
= 2 + 



12 

57i 2 

~6 



12 - 10A- 



h 2 



9i,j 



[I3q itj - 2g i _ij + gi_ 2 ,j 



which now needs three backward points to start. 

The appearance of the third r derivative of q"' in Eq. ( 10 ), which is strictly infinite at the 



step point, does not cause practical problems. Although not necessary, one can always 
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assume a smoothing of the potential at the cell boundary a la Becke, [28] reducing at 
the same time the mesh h, so that the error at that particular step point is negligible. 

In this way, at the cost of a bigger error 0(h 5 ) compared to the original Numerov 
formula and the introduction of a further backward point (three points i, i — 1 and 
% — 2 are now involved in (11 )), the three-dimensional discretized equation can be solved 



along the radial direction for all angles in an onion-like way, provided the expansion 
(|2| is performed at each new radial mesh point to calculate q it j. We use a log-linear 
mesh p = ar + (3\nr, to reduce numerical errors around the origin and the bounding 
sphere. [29J 

2.2. Matrix Numerov method 

It is well known that errors in the Numerov difference equation originating from the 
unidimensional differential equation 



d 2 



+ E 



1(1 + 1) 



- V(r) 



Pi(r) = 



(12) 



dr 2 r' 

grows exponentially when E — 1(1 + l)/r 2 — V(r) < 0. Therefore near the origin and in 
general for large r meshes and/or high I values the method is not suitable. This is also 
true for Eq. (jlj. To avoid this problem we use the so called Gaussian elimination for 
the difference equation [301 EH E2] • We notice that in the MT sphere lying inside the 
cell the AM expansion of the potential is regular and in general only few multipoles are 
appreciable. Therefore, by projecting onto Y L (r) we can rewrite Eq. ([I]) as [33] 



<h + l ^ r ^-E)X LL ,(r)+ I dvY L (v)V(v)P LI (v) 



dr 







(13) 



i.e. 



E 



L" 

= F(r)X(r) 



d 2 1(1 + 1) 
dr 2 r 2 



E 5 LL „ + V LL „(r) 



X VL »(r) 







where XiL'(r) = r Rlv(t), X is its transposed, 

V lv (t) = V VL (r) = [ drY L (r)V(r)Y L ,(r) 



(14) 



(15) 



and 



(F(r)) 



LL> 



F(r)) UL 
d 2 



+ l ^A-E) l ) Li/ V LL ,ir) 



ar z r z 

Eq. ( 14 ) is a system of coupled radial Schrodinger equations in matrix form that 



can be solved simultaneously for all L, L' components with appropriate initial conditions. 

The Numerov recursion relation for the matrix SE [33J is (notice the change of sign 
of the coefficient B for sake of later convenience) 



A,; + iXj +1 + BjXj + Ai_i~Ki_i — 



(17) 
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h 2 



Aj = 1 Pi 

12 



B, 



i ILL' 



5h* 
2+ Pi = 

6 1 

= v LI/ (n) + 



12 - 10A 

1(1 + 1) 



-E 



$l v 



where i is the generic point of the radial mesh. Its explicit matrix form is, 



/ 



A Bi A 2 

Ai B 2 A 3 



o \ 



\ 



o 



/X \ 
Xa 



V x M+ i / 







(19) 



Aj\/_i B M A M+ i J 

Since the regular solution has the boundary condition, X = we can rewrite this latter 
equation as 



Bi A 2 

Ai B 2 A 3 



o \ 



v 



o 



/Xi \ / 

x 2 



Am-i Bjv/ J 



\Xm J 








(20) 



\ -Aa/+iXm+i / 



This set of equations can be solved by performing forward Gaussian elimination near 
the origin, [301 



Di A 2 

D 2 A 3 



0\ 



O 



Dm-i Aa/ 



/Xi \ / 

x 2 



Xjif-i 
\X M / 










\ 



(21) 



V — AM+iXjvf+i / 



with 



Dx = B 1; 

D 2 = B 2 - AiDj-^a, 



D, = B, - A^Cv^A,, (i = l,...,M) (22) 

constituting a set of forward recurrence relations for the quantities Dj. In terms of these 
latter we finally obtain the following recurrence relations: 

Xj = -D^Ai+iXi+i, (i = 1, M) (23) 

the solution of which can be calculated backward starting from Xjv/+i = I, modulo a 
constant normalization matrix. As will be clear from the following, this initial matrix 
in practice will not be needed. Summarizing, our strategy to generate the cell basis 
functions Pl{ v ) — is the following. In a spherical domain around the origin, 
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inside which there are no discontinuities of the potential, we use the matrix Numerov 
method with Gaussian elimination (GE), since we can expand the potential in a well 
behaved series of spherical harmonics. We use the GE method to avoid the well 
know instability of the Numerov recursion relation near the origin when the angular 
momentum 1 is high, as mentioned above. As boundary conditions we use X = at the 
origin and Xm+i = I at the radius of the sphere, which is usually taken to coincide with 
the MT sphere inscribed in the cell. We then take the last three points of the solution 
so obtained to start the 3-d modified Numerov procedure outward across the potential 
discontinuity up to the cell bounding sphere. Since the local SE we are dealing with is 
an homogeneous equation, its solution is determined up to an arbitrary normalization 
constant (reflected in the second arbitrary condition Xm+i = I )• For the basis functions 
Pl{ y ) we never need such a constant, since only ratios of these functions appear in MS 
Theory, as clear in the following. Instead, when we compare with a definite solution, like 
in Fig. s [T] and [2] below, we need to provide the value of this solution at another point, 
usually the radius of the sphere. This means taking a value for Xm+i appropriate for 
this solution. 

It is also clear that the method can also be applied to generate by inward integration 
the irregular solutions needed to calculate the Green's function. 

This procedure is quite efficient and was tested against analytically solvable, 
separable model potentials, with and without shape truncation, obtaining very good 
results. In Ref. [22J we have shown the comparison between the analytical solution and 
the numerical one for certain directions in the special case of the truncated potential 
V(x, y,z) = a 0(\x\ — R c ) + bO(\y\ — R c ) + c9(\z\ — R c ), where 9 is the step function, R c 
= 3.78 au = 2.0 A and a = -0.05, b = -0.1, c = -0.15 Ryd, for an energy E = 0.3 
Ryd. For this comparison we used an Z max = 7 and a number of surface points on the 
Lebedev grid equal to 266. 

Fig. [I] shows the same comparison in the more stringent case of a discontinuity of 
the order of one Ryd. We took indeed a = —0.5, b = —1.0, c = —1.5 Ryd for E = 0.3 
Ryd. In this case, along the z direction, one can even observe a kink in the curvature 
of the solution, which is well reproduced and related to the discontinuity of the second 
derivative at the truncation value of R c = 2.0 A. As expected, for good agreement we 
had to increase / max up to 11 and take a number of Lebedev points equal to 1454. Notice 
here that the numerical method to generate the solution is really three-dimensional and 
does not take advantage of the separability of the analytical one. For this comparison we 
used the Matrix Numerov method with GE up to R c = 2.0 A = 3.78 au, then switched 
to 3-d Numerov. Due to the high potential step in this case, we took a number N of 
radial mesh points given by N = 834. 

In order to test the reliability of the method also in the case of a potential which is 
not truncated but varies substantially in sign and magnitude inside the defining region, 
we show in Fig. [2] the same comparison for the Mathieu functions, solution of the 
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Figure 1. Real and imaginary part of the numerical solution of the SE along z 
direction for the separable truncated potential given in the text, compared to the 
analytical one. (Color online) 



separable SE with periodic boundary conditions 



cP_ 

dx 2 dy 2 dz 2 



a 



y 



ip{x,y,z) 

- 2q x cos 2x + 2q y cos 2y + 2q z cos 2z)ip(x, y, z) 



(24) 



for the case where 

q x 

% 



1.0, a x = —0.455139, parity = even, period = ty 
0.3, a y = —0.044566, parity = even, period = n 
1.0, a z = +1.859110, parity = even, period = 2tt 

The energy eigenvalue is E = 1.359405, l max = 20 and the number of surface points is 
given by 1730. For the convenience of the reader the Mathieu functions are described 
in Appendix A The number of radial mesh points was equal to 250. Also in this case 



we employed both methods of integration with a switch radius of 2.45 au. 

In general the minimum number of surface points is chosen according to the rule 
that to integrate exactly J dQY L (Q) we need ~ (/ + l) 2 /3 points. [25J If we want to 
integrate the product of two Spherical Harmonics, I should be the sum of the individual 
l's; the same for a product three, etc... So for the 3-d Numerov we need to integrate 
only a product of two functions whose expansions are both truncated to a certain l max , 
therefore the number of points is (2/ max + l) 2 /3, whereas for the truncated separable 
potential and the Mathieu functions we have the product of three functions (one for each 
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Figure 2. Comparison of the particular Mathieu function described in the text with 
the one generated by matrix and 3-d Numerov method at four different angles. The 
switch radius was at 2.45 au. (Color online) 



space coordinate), so we need (3/ max + l) 2 /3 points. Moreover for the matrix Numerov 
method we again have (3Z max + l) 2 /3 due to the calculation of Vll>- 

2.3. A linear-logarithmic mesh 

To solve Eq. by Numerov procedure, there are several choices for the radial mesh. 
Due to the singularity of the potential near the origin we found that the best strategy 
in our case was to take a mixed logarithmic and linear mesh, as usual in atomic physics. 
[29j [301 E2] For non-MT calculation, especially with truncated potential, this mesh is 
the appropriate choice. In this case the new radial variable is 

p(r)=ar + (3 In r (25) 

with a and /3 constant. A constant mesh size of p can be taken in the interval 
Po < P < Pn- The initial value of po is chosen according to the empirical formula 

p = -0(10+ \nZ) (26) 
whereas the final p is defined as 

p N = po + Nh (27) 
so that a is given by 

a = (p N - (3\nr N )/r N (28) 
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taking j3, the mesh size h and the number of points N as input values. In the calculation 
of local basis functions we choose r N equal to the radius of the cell bounding sphere Rf,, 
(3 = 1.0 and put iV m 100-R&. Instead for Fig.s [T] and [2] we took respectively = 0.67 
and (3 = 0.05. The value of r = r(p) corresponding to a given value of p can be readily 
found by application of the Newton technique. [52] 

Following the change of variable in Eq. (25) the SE in Eq. (14) becomes 
F(p)Y(p) = where 



(F(p)) 



LL' 



dp 2 



+ [ a + 



P 



X 



1(1 + 1) + p {ar + f3/4){ar + py 



E) U 



>LL' 



+ [a + 



P 



V w (r) 



Y(p) 



P 



<\ f -X(r 



(29) 



where r = r(p). The same, mutatis mutandis, applies to Eq. (flj for the three- 
dimensional Numerov method. 

A comment is in order at this point. Strictly speaking by changing to the log- 
linear mesh the Eq.s (18 23) are not valid anymore, since r = cannot be realized (it 
would correspond to p = — oo) and therefore not explicitly implemented as boundary 
condition. By working out again the Gaussian elimination process when Xq is not zero, 



one arrives at the same equations (22), except that in the rhs term the zero of the i-th 



row is replaced by D { \ A Y , where Y Q is the value of (29) calculated at the first 
point pq. Now 



Yll'(p) 



'a + 



JL 

r{p) 



X LL ,(r(p)) = y/ar(p)+PyftUjR w (T{p)) (30) 



Since at the origin Rll> is diagonal in 1 and behaves like a spherical Bessel function, 
■\/r(p)R LL (r(p)) is of the order of 10~ 3 for 1=0 at p ~ 10~ 5 , 10 -8 for 1=1 etc.. We can 
therefore take Y = and use the simplified Gaussian elimination formulas Eq.s (18 23). 
We checked that this is a good approximation also for 1=0. 



3. Multiple scattering method for scattering and bound states 

3.1. Scattering states 



We begin by presenting the derivation of MSE for scattering states. In this case we seek 
a solution of the SE continuous in the whole space with its first derivatives, satisfying 
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the asymptotic boundary condition 



V»(r;k) 



k 



e ik - r + /(r;k)- 



ikr 



(31) 



167T 3 

where k is the photo-electron wave-vector and /(r; k) is the scattering amplitude. The 
factor (k/(16n 3 )) 2 takes into account the normalization of the scattering states to one 
state per Ryd. In the spirit of MST we partition the space in terms of non overlapping 
space-filling cells Qj with surfaces Sj and centers at Rj. Accordingly we partition the 
overall space potential V(r) into cell potentials, such that V(r) = Y^j v j{ r j)> where 
Vj(rj) takes the value of V(r) for r inside cell j and vanishes elsewhere. As clear from 
the following the zero value of the potential outside the cell is not necessary and can be 
replaced by any constant. The results will not depend on this particular value. Here and 
in the following i\, = r — Hj. The partition is assumed to satisfy the requirement that 
the shortest inter-cell vector Ry = R; — Rj joining the origins of the nearest neighbors 
cells % and j , is larger than any intra-cell vector r j or rj , when r is inside cell % or cell j . 
If necessary, empty cells can be introduced to satisfy this requirement. We also assume 
that there exists a finite neighborhood around the origin of each cell lying in the domain 
of the cell. [12] We then start from the following identity involving surface integrals in 
dr = da 

N 

£ / [G + (r' - r; k)V^(t; k) - ^(r; k)VG+(r' - r; «)] ■ n, do* 

,-=1 JSi 

= [ [G+(r' - r; K ) V^(r; k) - ^(r; k) VG + (r' - r; .)] ■ n d(T . (32) 

Here Q Q = w hh surface S D , centered at the origin o and Gq(t' — r; k) is the 

free Green's function with outgoing wave boundary conditions satisfying the equation 
(V 2 + k 2 ) Gq(t' — r; k) = 5(r' — r), where k 2 = E — Vq and Vq is an arbitrary constant 
equal to the assumed value of the cell potential outside the cell domain. The identity 



(32 )is valid for all r' lying in the neighborhood of the origin of each cell, since in this 
case the integrands are continuous with their first derivatives. In this context we shall 
use two distinct /c-vectors, defined respectively as k = \[E and k = \/E — Vq. This 
letter will appear in the expansion of the Green's function Gq(t' — r; k) by spherical 
functions. [16] Obviously k = k for Vq = 0. 



Equation (32), with the choice Vq = 0, can also be derived from the Lippmann- 
Schwinger equation 

^(r';k) =e ik - r ' + J Gq{y' - r; k) V(r) ^(r; k) d 3 r (33) 



satisfied by the scattering state (see Appendix C ). However we prefer to start from 



the identity Eq. (32) to take advantage of the arbitrariness of the constant Vq. For 



convenience of the reader we recall the expansions [13 



e ikr 



47r]T^y L (k)J L (r;£0 (34) 



L 
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G + (r' 



r;«) 



\ g«|r'-r| 



An \r' — r| 



Ji(r^; K)H£{Yi] k) (r- < r 4 

L 



(35) 
(36) 



Notice for future reference that in the = the solid spherical harmonics 

Ji(r; k) and H^(r; k) are to be understood as Ji(r) = r z Yj / (r)/(2/ + 1) and Hl(y) = 



-i-i 



y^(r), due to the well known expansion 

1 =Y — ^Y L (r)Yi<r) 



(37) 



which is the k — > limit of Eq.s (35) and (36) 



The heart of MST is the introduction of the functions $l(i"j; k) which inside cell j 
are local solutions of the SE with potential Vj(rj) behaving as Jl{ y j] k) for r 3 - —> 0. They 
form a complete set of basis functions such that the global scattering wave function can 
be locally expanded as [12] 

#r i; k) = 5>l(k)M*j;*) (38) 

L 

where we have underlined the k dependence of ^(r.,-; k) through its behavior at the 
origin. 

In order to find the asymptotic behavior in the outer region CQ we introduce the 
scattering functions in response to an exciting wave of angular momentum L: 

M*o, k) = J L (r o] k) + J G+(r - r' ; k) V(r' ) k) dV Q (39) 

Then, under the assumption of short range potentials (i.e. potentials that behave like 
l/r 1+e with positive e at great distances), letting r G — > oo and using expansion Eq. (36) 

Jl{*o] k) 

+ £fl+(r o ;A0 I J L ,(r' ;k)V(r> )MCk)d 3 r> 

J2 A °M \jl(to\ k) + £fl£(r ; k)T° L , 



in Eq. (39) we find 

^(r ;k)= Ys A °M 



(40) 



(41) 



where, in order to impose the asymptotic behavior in Eq. (31), A° L = i z Y^(k) (k/n) 1 ^ 2 
and T£ L , is the T-matrix for the whole cluster, equal to 



T° VL = J Jv{v - k) V(r ) ^ L (r G ; k) d 3 r 



(42) 
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In general for short range potentials decaying slowly, the asymptotic behavior in Eq. 



(41) is reached only at great distance from the origin of the coordinates (usually at the 



center of the atomic cluster under study). In order to limit the number of cells, so that 
the surface S Q just surrounds the cluster, we introduce the local solution 

$ L (r ;k) = J2R°L>L(ro)Y L ,(r ) (43) 
v 

in the outer region CQ Q , which can be obtained by inward integration of the SE starting 
from the appropriate asymptotic value H^,(r ; k). Therefore we take here 



ij(r o] k) = J2 \A° L (k)J L (r o] k) + $ L (r ; k)A° L (k) 



(44) 



Notice that the function <3?L(r ; k) in Eq. (43) (and consequently R° L , L {r ) ) is complex, 
unlike the functions $i,(rj;fc) that can be taken real, if the potential is real. If the 
potential has a Coulomb tail, the spherical Bessel and Hankel functions should be 
replaced by the corresponding regular and irregular solutions Fl{y ; k) and Gi(r ; k) of 
the radial SE with a Coulomb potential. Due to the possibility that the optical potential 
used for calculating the spectroscopic response functions be complex, it should be clear 
from the context that the formalism works also for complex energies and/or potentials. 
The extension to complex energies will come very handy when exploiting the analytic 
properties of the Green function. 



Insertion of the expressions Eq.s (38) and (44) into the identity Eq. (32) provides a 



set of algebraic equations (known as MSE) that determine the expansion coefficients 
A^(k) and the A° L (k) in such a way that the local representations are smoothly 
continuous across the common boundary of contiguous cells. Indeed, taking r' in the 
neighborhood of the origin of cell i ^ o, using the expansion Eq. (35) (since r is 



confined to lie on the cell surfaces), and putting to zero the coefficients of Jl(i^; K ) due 
to their linear independence, we readily arrive at the MST compatibility equations for 
the amplitudes A^(k) and A° L ,(k) 



jL' 



LL' 



,4,(k) = WTvAiw + Nr L ,A° L ,(k) 



(45) 



where 



TTlO 

n LV 



[#+(r 4 ; K)V^(r i; k) - $^(r 3 -; fc)Vfl£(r i; 



• n j d °j 



'Hi( r il ^)VJL'( r o; k) - JL'(r ; k)VHt{ri] k) \ ■ n da 



So 



[H+( r i;K)V$L>(r ;k) -$i/(r ;fe)VflJ(r f ;«)] -n da 



So 



A further set equation is obtained by taking r' inside the outer region CQ Q , using 



the expansion Eq. (36) (remembering that r D < r^, since r Q lies on S ). By putting to 

(46) 



zero the coefficients of H^(r' ; k) we obtain 



J2K J L'AUk) = J2 [MlUUk) + N? L ,A° L ,{k) 

jL' V 
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where 



K lv = / [ J l{*o\ K)V$L'(r j; k) - <$>L>{rj] k)VJ L (v ; k) ] ■ n, daj 
Js 3 

M° L ° L , = 5 LL , / [J L (r ; K)VJ L ,(r ; k) - J L ,(r - k)VJ L (r o] k) ] • n da 
J s 

N° L °v= / [J L (r ;K)V^ L >(r ;k) -$ F (r ;fc)VJ L (r ; K )] -n da 
Js„ 



From the above derivation it is clear that the set of equations in Eq.s (45) and (46) 
determines the amplitudes A^(k) and A° L {\s) independently of the constant Vo, since 
the identity Eq. (32) is valid whatever Vq. In general this will be true only if the L- 
expansion is not truncated, whereas there will be a more or less pronounced dependence 
according to the degree of convergence of the truncated expansion. In general, the lesser 
the potential jump at the boundaries of the various cells the faster the convergence. 

Notice that these equations remain valid, with no restriction on the sums over L, 
even in the case k — 0, provided Jl and H L are replaced by J L and H L , due to the 
expansion Eq. (37) of the zero energy limit of the free Green's Function. 

The usual derivation of the MSE now proceeds by re-expanding H^^^k) and 
Jl{y ]k) around center j under the geometrical conditions stated at the beginning of 
this section, by use of the equations [TBI EE] 

k) = Y, &LM*j\ «) (RiJ > rj) (47) 

J L (r ; k) = ) j J^ L ,J L >(rj] k) (no cond.) (48) 
v 

H+(r i; = JwHUto; k) (r > R lo ) (49) 
v 

where G l [ L , are the free electron propagator in the site and angular momentum basis ( 
KKR real space structure factors) given by 

G% =^J2 C ^ L '> L ") k) (50) 

L" 

and J % l L , is the translation operator 

4v = ^J2 C ( L > L '> L ") J v( R v> *) (51) 

L" 

In these formulas the quantities C(L, L'; L") are the real basis Gaunt coefficients given 
by 

C(L,L';L") = J Y L {tt)Y L ,(tt)Y L ,,(tt)dtt (52) 

In the following we shall also need the quantity 

Nlv =^J2C(L,L';L")z l ~ l ' +l " N L/ ,(R tJ ; K ) (53) 
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Unfortunately the re-expansions Eq.s (47), (48) and (49) introduce further 



expansion parameters into the theory (with related convergence problems) that are 
actually unnecessary, as shown below. 

We in fact observe that the integrals over the surfaces of the various cells j can 
be calculated over the surfaces of the corresponding bounding spheres (with radius RQ 
by application of the Green's theorem, since both H^(r;K) and $L(r;k) satisfy the 
Helmholtz equation (V 2 + k 2 ) F(r) = outside the domain of the cell. We then use the 
following relations 



Y L ,(r 3 )VH L (r<) ■ n, da, = (Rjf G% -ji^K) 

dR J b 



(54) 



(55) 



which are exact for all L provided |r^ — r,-| = Rij > r,- for r lying on the surface Sj. This 



is a consequence of the fact that under this condition the series in Eq. (47) converges 



absolutely and uniformly in the entire angular domain, as shown in |Appendix B[ Eq 



(B.7) and can therefore be integrated term by term. This property is also true for 



the series derived with respect to r. Even though not necessary, we also checked the 



numerical equality of both sides of Eq.s (54) and (55) for various values of L, V 



in 



Similarly, since the series in Eq. (|49|) converges uniformly and absolutely, as shown 

(56) 



Appendix B Eq. ( B.8[ ), we also find 

Y L ,{r ) Ht(r i; K ) da, = {R° b f J L ° L , h+(nR° b ) 



So 



Y L ,(i )V H+in) • n, da, = {R° b f T L ° L , — h+( K R° b 



(57) 

'S„ U£X b 

provided Ri Q < R b , where R b is the bounding sphere of the outer region CQ . Therefore 
R b should be bigger than any R io . 



Finally, due to the absolute and uniform convergence of the series in Eq. (48) 
without conditions, we find the following relations 



M?,-) J L (r,; k) da, = (i^) 2 J% J^kR, 



YjAt^V JAr l} 



n,- da,- 



{R j b ?J ij 



dR J b 



(58) 



(59) 



By inserting in Eq. (45 ) the expression for the basis functions expanded in spherical 



harmonics (we shall suppress the site indices whenever a relation refers to both sites i 
and site o) 



* L (T;k) = y %2R L , L (r)Y L ,(r) 



(60) 
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remembering that this expansion is uniformly convergent in the angular domain [24J and 



using the relations Eq.s (54)- (59) we finally obtain, under the partitioning conditions 



specified at the beginning of Section [3j 

££i L ,Ai,(k) + ]r 

V j,L',L 



^LL"^L"L'^L'(^) 



LL> 



M£^,(k) + ^££ a ,,^„(k) 

L" 



(61) 



where we have put E° LtLII = Nl° L „, the quantities M™ L and NZ? L „ being the same as 



those following Eq.(45), calculated with replaced by r . 
Similarly, putting S° LL , = Nf^i, for Eq. ( |46| we find 

j,L',L" 

In the above equations we have defined the quantities 



\M° L ° L ,Al,(k)5 LL , + S° LL ,A L ,(k) 



(62) 

(63) 
(64) 

for the cells flj and for the outer region CQ Q . The Wronskians W[f, g] = fg' — gf are 
calculated at Rl and Rl respectively and reduce to diagonal matrices for MT potentials. 



Ell> 
Sll' 



(Ri 



?W[-iKht,R LL , 
) 2 W\j h R w \ 



Equations (61 ) and (62) look formally similar to the usual MSE. However we notice 



that due to the relations Eq.s (54)-(59) there are only two expansion parameters in the 



theory. They are related to the AM components of Rl'l in the expansion Eq. (60) 
in cell j and in the outer region CQ Q . No convergence constraints related to the re- 
expansion of the various spherical Bessel and Hankel functions around a different origin 



Eq.s (47)-(49) are present. 



It is interesting to note that the truncation value for both indices is the same 
and corresponds to the classical relation Z max = kR J b , where R{ is the radius of the 
bounding sphere of the cell at site j. This is true for the index L, which reminds that 
the basis function $l is normalized like ji(kr)YL near the origin. Due to the properties 
of the spherical Bessel functions, when I ^> kR J b , $^ becomes very small inside the cell, 



decreasing like [(2/ + l)!!] -1 . Therefore his weight in the expansion Eq. (60) will be 
negligible. The other index L', as will be clear from the following, measures the response 
of the truncated potential inside the cell to an incident wave Jy of angular momentum 
L' . Due to the same argument as above, familiar to scattering theory, the scattering 

for 



see Eq. (B.10) in Appendix B 



matrix T[, L will decrease like [(2/ + 1)!!(2/' + 1)!!]" 
/, /' 3> kR 3 b ). As a consequence and S" J can be considered square matrices. In the case 
of the outer sphere region CQ , the situation is inverted, the index L being related to 
the response of the entire cluster to an incident wave of angular momentum L, whereas 
the index L' corresponds to the number of AM waves mixed in by the potential not 
only inside Q but also in CQ Q . The two indices have the same truncation / max = kR b , 
provided we take R b as the radius of the sphere that contains the region of space where 
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the potential is substantially different from zero. This conclusion is reinforced by the 
observation that one can cover this same region by empty cells. 

Up to this point we have assumed that Vq ^ and derived consequently the MSE, 
having in mind the possibility to check the rate of convergence of the L-expansion. 
However in the continuum case one usually works under the assumption that Vq = 0. 
In this case the Eq.s (61) and (62) simplify considerably in the case of short range 
potentials. Since now k = K, we use the relation 

[H£{r ; k)VJ L {r - k) - J L (r ; k)VH+{r ; k) } • n, da = -8 LU (65) 

So 



(66) 



so that in Eq. (61) Ml° L = —1, and in Eq. (62) M£^ = 0. Moreover one easily finds 
that 

£iMk) Jtv = i l Y L Qz) ^°Jl = I[(k) 
v 

which is obtained from Eq. ( |51~j ) by observing that 

J2C(L,L';L")Y L ,{n) = Y L (Q)Y L „(n) 

L> 

Then the two sets of equations assume the simpler form 
J>i L Ml,(k) + ]T G i i L „S J L „ L ,A j L ,(k) 

V j,L',L" 

^Y, Jl ^ E l'L''A° L!! (k) = -PM 



L'L'' 



3+o 



(67) 
(68) 



J^Si f/LI A j Lf (k)-J2s°LL'A o LI (k) = 

j,L',L" L< 

The fact that E and 5* can be taken to be square matrices leads to another interesting 
form of the MSE. Under the assumption that Det S ^ 0, we can introduce new 
amplitudes 

B L (k) = J2SLL'A L ,(k) (69) 
v 

which is equivalent to using new basis functions $l related to $l by the relation 



where 5* is the transposed of the matrix S. 
Defining the quantities 

(T*) -1 = -E^S 1 )- 1 
T° = -E^S )- 1 



(71) 
(72) 
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(notice the asymmetry between sites i and site 6) we can write Eqs. (67) and (68) as 

^(r)-i,si,(k) -f>^,i?£,(k) 

^J^Tla»^,(k) = /i(k) (73) 



L'L' 1 



3+o 



(74) 



^«(k)-B»(k) = 

j,L' 

The meaning of the amplitudes -Bi(k) is immediately found from these equations if we 
consider only a single truncated potential at center i. In this case T° = 0, since now the 
asymptotic behavior is given by Eq. Mil), and B L (k) = A L (k) = T£ L ,A L , where T£ L , 



is the T-matrix of the potential. Therefore Eqs. ( 73 ) and ( 74 ) tell us that T l LL , = T 



LV 



As a consequence B 1 l (\l) is the scattering amplitude of angular momentum L in response 
to an exciting plane wave of wave vector k. Moreover, we find that T l = — S^i*?) -1 is 
symmetric in the AM indices (remember that we use a real spherical harmonics basis), 
a fact already known from general scattering theory. This is a consequence of the fact 
that SE^ 1 is a symmetric matrix. [31] 

In the case of many cells, it is expedient to work only in terms of the cell amplitudes 



B l L ,(k). Inserting into Eq. (73) the expression for I?£,(k) given by Eq. (74) we obtain 

^(r) L i,5i,(k) -£G^i,(k) 



j,L' 



I'M (75) 

jL' AA' 

Introducing r, the inverse of the multiple scattering matrix M = T _1 — G — JT° J 

(T- 1 -G- JTJ)- 1 (76) 



T 



known as the scattering path operator [13], we derive from Eq. (75) that 



51(k) = E^ J i'( k ) 



(77) 



jL' 



If we insert this expression in Eq. (74) and remember that by definition -B£(k) ; 
Y2l' ^■ll'^ l l'i we eas ily find for the cluster T-matrix 

J LA T AA' Ja'L' ( 7£ 



I LL' 



ij AA' 



Since the matrices G and J are also symmetric (see definitions Eq.s (50) and (51)), we 
find that r is likewise symmetric, implying the symmetry of T£, L , again in keeping with 
scattering theory. This quantity represents indeed for the whole cluster the scattering 
amplitude into a spherical wave of angular momentum L in response to an exciting wave 
of AM L' and is needed for example in electron molecular scattering. [33j Finally Eq. 
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(77) shows that the quantities B l L (k) are scattering amplitudes for the cluster, for which 



the generalized optical theorem holds (for real potentials) 

1 



16] (see Appendix D ) 



dksi(k) [Bi,(k)j 



^ r 



7T 



LL' 



(79) 



This relation is very important, since it establishes the connection between the photo- 
emission and the photo- absorption cross section, as shown in Appendix E As it will 
turn out, — Sr|^ is proportional to the L-projected density of states onto site i. 

In the case of one single cell located at site i, by construction the solutions inside 
and outside the cell are continuously smooth so that, remembering that by definition 



Tl L , = T° T/ , for r. 



LL'-> 



R l b we have, neglecting for simplicity from now on the k 



dependence of the local solutions, 

^i?l(k)$ L (r,) = ^i2(k) 



J L (r ;k) + J2HUr ;k)T L/1 



(80) 



L L 

Using Eq. ( |77| for a single site and equating the coefficients of A° L (k) we find at the 
bounding sphere the relation 



J L (r ;k) + J2HU^;k)T L/1 



(81) 



implying that the basis functions $ L are scattering functions, obeying the Lippmann- 
Schwinger equation for the cell potential. Therefore, introducing new expansion 
coefficients Cx(k) such that locally 

^(r;k) = ^C L (k)$ L (r) (82) 



and repeating the steps leading to the MSE in this new basis, we obtain 

c[(k) - J2 G^T|, /L ,ci,(k)-^jr , L'^(k) = /i(k) 

3, L'L" V 

E ^i" T l»L^i'(k) = E( r )^^'( k ) 

3, L'L" V 



(83) 



34) 



Comparing these equations with the previous ones in Eqs. (73), (74) and (67), (68) we 
immediately find the relations 



L> 

C L (k) = J2 E LL'A L ,{k) 



(86) 
(87) 
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In the present approach, the three forms of pair of equations (67)-(68), (73)-(74) and 
(83)- (84) are equivalent and lead to the same result. 

The pair of equations (83)- (84) are quite important, since they provide the formal 
justification that in MST one can work with square matrices, provided that the only 
indexes appearing in the theory are those of the radial functions Riy{r). This is a 
consequence of the relation (B.9) of Appendix B (second equation) and the fact that 
the matrix elements T^y have a common truncation parameter Z max . In fact, since 
Tr(T^T) < oo due to the asymptotic behavior of the Tlv matrix elements given by Eq. 
(B.10) in the same Appendix, one can safely define an inverse for the matrix Tll'(E) 
(except at poles on the negative energy axis) and pass from one representation to the 
other. In particular one can pass from the set (83)- (84) to the set (67)- (68). In the 
traditional derivation of MS equations, that does not rely on the relations (54)-(59) but 
hinges on the re-expansion formulas (47)-(49), this equivalence does not hold. In fact the 
need to saturate the "internal" sum over L" coming from the re-expansion introduces 
a further expansion parameter and therefore rectangular matrices into the theory. This 
feature makes it impossible to define a T-matrix and to write a closed form for the 
GF, loosing all the advantages of MST over other methods. This drawback has been 
avoided in our approach, since in each step of the derivation of the MS equations we 
have shown that the introduction of summation indices other than those present in the 
radial functions -Rl'l(^) is unnecessary. 

Another useful consequence of the fact that the theory can be cast in terms of 
square matrices is the possibility to exploit the point symmetry of the cluster under 
study. Even though many authors have treated the problem of how to symmetrize the 
MSE, this was done in the framework of the MT theory, where the cell T-matrices are 
diagonal in the AM indices. New features appear in the more general case (in particular 
how to calculate the symmetrized version of the T^y matrices) and Appendix F deals 
with this situation. Needless to say, we checked in all applications that the symmetrized 
and unsymmetrized version of the theory gave the same results. The application of 
the symmetrization procedure to Green's Functions or to periodic systems is rather 
straightforward. 

As already anticipated in the introduction, one of the major advantages of MST is 
the direct access to the Green's Function of the system. Having explicit expressions for 
this quantity is of the utmost importance both for writing down spectroscopic response 
functions (see Ref. [35]) and for the calculation of ground state properties through 
contour integration in the complex energy plane (see e.g. Ref. [9] and references therein). 

The GF is solution of the Schrodinger equation with a source term 

(V 2 + E - V (r) ) G(r, r'; E) = 5(r - r'). (88) 

In the framework of MST and for general (possibly complex) potentials, the solution 
of this equation in the case of a finite cluster can be written as [131 136] 



G(r i ,i> i ;E)= ($(r< 



6ij T 



+ ^•(<l>(r < )|T'|^(r' > )) 



S9) 
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where r< (r>) indicates the lesser (the greater) between and r[. The function \I/(r) is 
the irregular solution in cell i that matches smoothly to H^,(r) at R l b . For short we have 
saturated the sum over the angular momentum indices using a bra and ket notation (e. 
9-) 

(5(r 4 )|T*|$(r})> =X;^( r *)^/^W) (90) 

LL> 

Moreover, for simplicity of presentation we have assumed no contribution from the outer 
region potential (i.e. T° = 0) allowing empty cells to cover the volume Q Q up to the 



point at which the asymptotic behavior in Eq. (41 ) starts to be valid. The modifications 
needed in the case T° ^ are obvious. In the case of a crystal we have to work in Fourier 
space [9]. 

Now, from Eq. (l8T|) written as 



$ L (r.) = -JL>(r ; k){T- l y L , L + tf+(r D ; k) (91) 



v 



by continuity we derive inside cell i the relation 

ttfo) = A L'(r 4 ; k)(T~% L + V L (r f , k) (92) 
v 

where Ajy(rj) is the irregular function joining smoothly to Jl'(t ; k) at R l b . Therefore 
the Green's function takes the form 

G(r l ,r' 1 ;E)= ($(r,) | \ $(r$) ) 

- ^•($(r < )|A(r' > )) (93) 

For real potentials, both and are real, so that the singular atomic term 
does not contribute to the imaginary part of the GF. In this case the quantity 
— S J n G(r, r; E) d 3 r = — 3 t™ l (E) J n $ L (r) d 3 r is the projected density of states 
on site i at energy E, expressed as a sum of the partial densities of type L. This relation 
(not r-integrated) constitutes the basis for calculating the system density by contour 
integration in the complex energy plane. 

Alternative forms of the GF that are independent of the normalization of the local 
solutions $z,( r i) can be easily obtained in terms of the S and E matrix. For example 
we have 

Giv^-E) 

= -<*(r«)| {{[~SE+~SGS\- 1 f-8 ij {[~SEr 1 f}\H^)) 

- 5 l3 (<$>(v < )\E- l \^(r' > )) (94) 

which is seen to reduce to the following expression, remembering the definition of | $ ), 

G(r u r' J ;E)= mn)\([I- GT]" 1 G) ij |$(r})> 

- <%<t(r<)|*(r'>)) (95) 
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Indeed from the relation, 

(A + B)- 1 - A' 1 = (A + BY 1 (A - (A + B)) A~ l 
= -{A + B)~ 1 BA- 1 
= - {B~ 1 A + l)-^- 1 

= - (AB^A + A)- 1 (96) 

we find 

[SE + SGS]- 1 -[SE]- 1 = -[SE+ SEISGS^SE}- 1 

= - [SE+ S EIGS]- 1 E]- 1 
= - E- 1 [3 + SEIGS]- 1 ]- 1 
= - E~ x [ S + ESIGS]- 1 ]- 1 
= E~ l [T -G- 1 ]- 1 ]}- 1 
= — E~ x [I — GT]~ l G E~ l (97) 

taking into account that S E = E S and T = —SE^ 1 = —E~ 1 S. All these forms are 
equivalent as long as we can treat the matrices S and E as square. 



3.2. Bound states 



Even though the essential of this section has been presented in a conference proceedings 
[22] . we feel that for the sake of completeness of presentation and convenience of the 
reader it should be repeated here. 

The MSE in the case of bound states can be derived from those for scattering 



states, by simply eliminating the exciting plane wave in Eq. (109) and taking the 
analytical continuation to negative energies in free Green's function Gq(t' — r;k), in 
order to impose the boundary condition of decaying waves when r' — > oo. In this case 
the Lippmann-Schwinger equation reduces to the eigenvalue equation 



Gj(r'-r;fc)^(r)^(r) d 3 r 



(98) 



where we have dropped the label k in the wave function ip(r'). Since the expansion 



of Gq(t' — r;k) in terms of spherical Bessel and Hankel functions in Eqs. (35) 



and (36) remain valid under the analytical continuation to negative energies, so that 
k = \/E = i\/\E\ = ry, we see that ip(r') behaves like e lkr ' jr' = e~ ir ' /r' for r' — > oo. 
We remind that 



ht(kr) = -r l KlM;K(kr) = -r\-iy Kf( ir ) 



ji(kr) = i l Ii^yr); ni(kr) = i 



(99) 



where // is the modified Bessel and Kf, Kf the modified Hankel functions of first and 



second kind, respectively. Not only the expansions in Eqs. (35) and (36), but also 



the re-expansion relations in Eqs. (47), (48) and (49) remain valid under analytical 



continuation with the same convergence properties (see Appendix B). This fact implies 
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that we can derive the MSE for bound states following the same patterns as for scattering 
states, except that now the behavior of the wave function in the outer region CQ is 

#r )= $>£<l£(r ) 

L 

= J2 A °Lj2 R °L'L(ro)Y L ,(i ) (100) 

L V 

The functions $2( r o) are now rea l an d can easily be found by inward integration in the 
outer region starting from an asymptotic WKB solution properly normalized, e.g. like 
[(2/ + l)!!] -1 . 

Working with the Bl amplitudes we easily arrive at the following condition for the 
existence of a bound state 

E I - (1 - Si,) G LL> ~ E ji LL^VV'4'V \ BL = (101) 

jL' I L'L" ) 

which is the same as Eq. (75), except that the exciting plane wave term I^(k) and the 
k dependence have been dropped. Notice that we have kept the arbitrariness of Vq in 
the free Green's function, in order to check that the eigenvalues do not depend on it. In 
the spirit of the analytical continuation, we have a definite rule on how to calculate the 
various quantities function of K. 
We now define 

C LU = {R b ) 2 W[n h R LU ] (102) 



so that, remembering Eq. (72) 

K -i( T i)-i = ( K jyi + i = _ C j( S jyi + i ( 103 ) 

k~ x T° = K° + i = -C^S )- 1 +i (104) 

Moreover we observe that 

«r x G% = N^ L , - iJ*, (105) 



where Nfj L , is defined in Eq. (53) and that ^2 L „ Jll"^l"L' = J % ll'i s i nce J is the 
translational operator. Substituting these relations into Eq. (101) and eliminating the 
common factor k _1 we finally find 

E j - C 1 - *v) N ll' - E J LL'K>v>J°Li> \Bh = (106) 

jL' L L'L" ) 

The generic (LL')-element of this MS matrix is either real for real k {E — Vq > 0) or 
proportional to i l ~ l ' +1 for imaginary k (E — Vo < 0). Indeed, due to the relations Eq. 
g, putting for short Ki = K} + Kf}/2, we easily find that 

N LL , =4n i l - l ' +1 J2C(L,L' ] L'')(-l) 1 '' ^^IkIR^Yl^R^) 



L" 



kIl' = -i^msr 1 ]^ 
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where C_ and S_ are denned in terms of the modified spherical Bessel and Neumann 
functions as the corresponding quantities. 

Therefore the condition for a bound state becomes Det M = 0, where M is the MS 



matrix in Eq. (106) after a unitary transformation that eliminates the imaginary 
factors. In the practical numerical implementation we find the zeros of the determinant 
of Det (KM), excluding the spurious solutions coming from the zeros of Det S. In this 



form, the procedure is equivalent to finding the poles of the GF in the form Eq. (95) 
on the real negative it should be. Still numerical instabilities might come from 

the inverse of S° present in the contribution of the outer sphere region. This unwanted 
feature could be eliminated by working with the A L , instead of the B L amplitudes. 
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Figure 3. Partitioning of the space for the hydrogen molecular ion with no empty 
cells. 

We applied the theory above to find the exact eigenvalues of the hydrogen molecular 
ion, since this test is considered rather stringent for the validity of the theory due to 
rapid variation of the potential in the molecular region and to the awkward geometry of 
the cells. In this case we partition the space in three regions, as illustrated in Fig. |3j two 
truncated spheres around the protons with a radius of 1.72 a.u. corresponding to cells 
Qi and flu and an external region labeled Qui, corresponding to the complementary 
domain CQ Q . The bounding sphere of this latter is represented by the dashed circle with 
radius 1.4 a.u., bigger than one half the distance of the protons, as discussed after Eq. 
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Table 1. Eigenvalues of the hydrogen molecular ion (in Ryd) 



Mol. 


n 1 m 


Exact 


Smith & 


Foulis [34] 


22 EC 


22 EC 


No EC 


No EC 


orb. 






Johnsonl37l 




V =-1.90 


V = 


V =-1.90 


V = 


lai g 


100 


-2.20525 


-2.0716 


-2.18973 


-2.20522 


-2.2055 


-2.2050 


-2.2048 




200 


-0.72173 


-0.70738 


-0.72093 


-0.723 


-0.724 


-0.731 


-0.726 


3ai g 


320 


-0.47155 


-0.45574 


-0.47102 


-0.4727 


-0.478 


-0.476 


-0.474 


4a lg 


300 


-0.35536 


-0.34859 


-0.35525 


-0.356 


-0.3550 


-0.357 


-0.356 




2 1 


-1.33507 


-1.2868 


-1.33426 


-1.3348 


-1.3348 


-1.3342 


-1.3343 


2a2 U 


3 1 


-0.51083 


-0.49722 


-0.51085 


-0.51072 


-0.5105 


-0.5104 


-0.5104 


3a2« 


4 1 


-0.27463 


-0.26979 


-0.27466 


-0.27469 


-0.2742 


-0.2745 


-0.2745 


4a 2 „ 


430 


-0.25329 


-0.24997 


-0.25329 


-0.254 


-0.2536 


-0.2541 


-0.25301 


le i 9 


32 1 


-0.45340 


-0.44646 


-0.45333 


-0.4545 


-0.45332 


-0.455 


-0.455 


lei„ 


2 1 1 


-0.85755 


-0.88866 


-0.85585 


-0.85754 


-0.8561 


-0.870 


-0.858 



(57). By calling the region outside this circle Cf2&, the potential is taken to be zero (or 
constant) into the intersection of this domain with cells Qi and Qu, and equal to the 
value of the true potential in the intersection with CQ Q . We also did a calculation with 
the two atomic cells, 22 empty cells surrounding them, plus an external region. 

It should be noticed that treatment of bound state is done here in analogy to the 
X-a MST method [3J, since we intend to put the theory to a severe test concerning 
the independence of the eigenvalues from the value of the interstitial constant Vq and 
the partitioning of the space. More modern techniques that avoid finding eigenvalues 
and eigenstates of the molecular cluster in the course of a SCF-iteration, exploit the 
analyticity of the GF through a contour integration in the complex energy plane to find 



directly the molecular density, as mentioned in the introduction and in Section 3.1 

Our findings are listed into Table ([!]) and compared with the exact results. The last 
two columns show the eigenvalues obtained with two different values of Vo, respectively 
equal to -1.90 Ryd and 0, showing the 'quasi' independence of the results from the 
constant interstitial value Vq. The columns with the label '22 EC refer to the calculation 
with two atomic cells, 22 empty cells and an external region, showing the 'quasi' 
independence of the result from the partitioning mode of the space. We attribute the 
slight dependence of the eigenvalues on V and the partitioning mode to the numerical 
instabilities mentioned above and the L-truncation of the matrices. 

The column labeled 'Smith & Johnson' refers to the calculation by Smith and 
Johnson [37] in the MT approximation, whereas the one labeled 'Foulis' quotes the 
result by Foulis [31] obtained within the distorted wave approximation. 
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4. Convergence of Full Potential Multiple Scattering Theory 

The inversion of the MS matrix becomes computationally heavy at high photoelectron 
energies because of the large number of angular momenta involved, since Z ma x ~ kRb- A 
common way to circumvent this difficulty is to invert the MS matrix by series, whereby 

(T- 1 - G)- 1 = T ^(GT) n . ( 107 ) 

n 

While this series is absolutely convergent for nonoverlapping MT spheres, provided the 
spectral radius of the matrix GT is less than one [38J, it is known to diverge for the case 



of space- filling cells. This is easily seen by using the inequality (B.ll) in Appendix B 
putting I = I' Z max , whereby 

/7Rh\ 21+1 

\GuT u \ « # 6 (^) (108) 

which signals the divergence of the matrix element (GT)u (for space-filling cells 2R^ > 
Rij, at least for nearest neighbors). 



However due to the behavior shown by Eq. (108) there is a widespread belief that 
the procedure of inverting exactly an / truncated MS matrix and then letting I go to oo 
does not converge in the case of space-filling cells. We shall show in the following that 
this is not so, provided a slight modification of the free propagator G is adopted. 

In order to illustrate our point, let us start by solving the Lippmann-Schwinger 
equation using the theory of the integral equations, before applying MST. 

V>(r';k) = e ikr ' + J Gq{t' - r; k) V(r) ip(r; k) d 3 r (109) 

We cannot use Fredholm theory, since the Kernel for this integral equation 

K(r',r) = -i-^^V(r) (110) 

is such that 

Tr(K f K) = / / drdr'K*(y,r)K(y,r) 




Air J J J \r' — r| 2 

and obviously diverges. 

However a solution for this problem can be found by following the argument of 
section 10.3, page 280 of Ref. 39 in the paper (Newton). We multiply the Lippmann- 



Schwinger equation Eq. (109) by iV^r')! 1 ' 2 and write V(r) = |V A (r)|w(r), where v(r) is 



a sign factor, equal to +1 where the potential is positive and to —1 where it is negative. 
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Then we obtain 

^(r';k)^|V(r')r /2 ^(r';k) 



\V(r')\ 



All/2 ik-r' 



+ \V(r')\ 1 ' 2 j G?3-(r'-r;ife)|y(r)| 1 / 2 t ,(r)V i (r;k)d 3 r (112) 



The Kernel for this integral equation is given by 



K fl (r',r) 



L\V(v')\ 1/2 - iTfvU' - 



An 



\V(r)\V 2 v(r) 



(113) 



whereby 



Tr(K+K s ) = j y'drdr / K*(r , ,r)K(r , ,r) 

jV(v)\V(v')\ 



< 



1 

All 

1 
47T 




drdr 




drdr 



,\V(r)\\V(r>)\ 



which is finite for a large class of potentials (including the molecular ones), so that the 
kernel K s is of the Hilbert-Schmidt type and Fredholm theorem for /^-kernels can be 



applied. Once the solution ^(r'jk) is found, we can obtain the solution of Eq. (109) 
simply by dividing it by ^(r')) 1 / 2 , except at points for which |K(r')| 1//2 = 0, where it 
can be defined by continuity 



Now, let us apply MST to Eq. (109) using the scattering wave functions $L(r) 



in Eq. (81) as local basis functions. We transform this Lippmann-Schwinger equation 



into a set of algebraic equations of infinite dimensions for the coefficients C^(k) in the 



expansion (82) 



ci(k) - Gi L ,,Ti, IL ,ci,{k)=ii{k) 

j,L'L" 



(115) 



where, in comparison with Eq. (83), for simplicity we have neglected the outer region 



C£l , which we can always assume to be covered by a set of empty cells. In matricial 
form we have, putting K = GT and calling A the term in the rhs, 

(I-K)C = A (116) 

The matrix K here is not an operator of the Hilbert-Schmidt type, since Tr (K^K) 



diverges, due to Eq. (108) and in keeping with Eq. (111). However, following the 



procedure used above in passing from Eq. (109) to (112) and introducing new vector 



components C = T l / 2 C with a new inhomogeneous term A' = T 1 / 2 'A, we transform 
this equation into a new one 

(I-K S )C' = A' (117) 

where 

K T l/2Q T l/2 , m) 



Full- Potential Multiple Scattering Theory with Space-Filling Cells for bound and continuum states31 



The square root of the matrix T is defined in the usual way, by first diagonalizing it 
with a similarity transformation S, taking the square root of the diagonal elements and 
then performing the same transformation on these letters. In formulas, if A = STS^ 1 , 
then T 1 / 2 = SA. 1 / 2 ^ -1 , so that T 1 / 2 T 1 / 2 = T. There is no danger in performing these 
operations with the infinite matrix T, since Tr (T^T) < oo, as can be seen from the 
asymptotic behavior of its matrix element in Appendix B Eq. (B.10). Hence the 
limiting procedures are well defined. 

By virtue of Eq. ( 114), Appendix G shows that the kernel K s here is of the Hilbert- 
Schmidt type (i.e. Tr (K\K S ) is finite). As is well known |39j this letter is the condition 
for the existence of the determinant \I — K s \ necessary to define its inverse, since by 
Hadamard's inequality, for any finite L max , one has 



\I-K, 



| 2 ^ J~j£max 



s LL> 



119) 



and in the limit L max — > oo the infinite product will converge if \{Ks)ll'\ 2 = 

Tr (K\K S ) < N < oo. jlO] 

This means that the process of truncating the matrix / — K s to a certain Z max 
and then taking the inverse, converges absolutely in the limit Z max — > oo. Once C is 
obtained, C = T~ 1 / 2 C, thus solving the original problem. Moreover the scattering path 
operator r (76) is given by 

T l/2 (/ 



T 



K^T 1 ' 2 

There is another way to solve Eq. 
writing 

(/ - K s y l = J2(KsT 

n 



(120) 



117, by expanding (I — K s ) in series, i.e. 



121) 



However, even if the kernel K s is of the Hilbert-Schmidt type but Tr (K\K S ) > 1, 
the series diverges, whereas the process of truncating and taking the inverse always 
converges. It goes without saying that the series ^2 n K n is always divergent, since 
Tr (K*K) is infinite. Therefore the series expansion procedure is not always a viable 
method to find the inverse of a matrix of the type (I - A) . 

In practical numerical applications one does not have to worry about modifying 
the structure constants according to Eq. (G.5) since, for the cell geometries ordinarily 
encountered in the applications (see the restrictions described at the beginning of 
Section [3j), /-convergence in the /-truncation procedure of the MS matrix shows up 



much earlier than what predicted by the onset of divergence in Eq. (G.4), written with 
the unmodified structure constants G^A'- We already found this out in GeCl 4 [2Tj . 
where in the first 20 eV an Z max = 3 was sufficient to reproduce all spectral features, 
which did not change by increasing / up to 10. Similar results were found for other 
compounds. 

In this context we did a more stringent test for the diatomic molecule formed 
by two inter-penetrating nonequivalent spheres with 40% overlap, with centers on the 
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z-axis. We calculated the K-edge ^-polarized cross section for the a state (m = 0) up 
to Z max = 60 in the energy range —4.0 ~ 20.0 eV, using the kernel K = GT and found 



a convergent behaviour for the lhs of Eq. (107) (see Fig. |4J). 

The fact that the full inversion of the MS matrix (/ — GT) is stable in this case 
up to Z max = 60 is clearly not of general validity, although indicative of the behavior 
of the theory. Going to higher values of Z max is not easy, because the lack of Lebedev 
integration formulas for a number of surface points > 6, 000 prevents us to access such 
values. Already the slight discrepancy of the Z max = 60 curve in Fig. [4] with the previous 
ones (barely visible) is the sign that ~ 6,000 Lebedev points are barely sufficient in this 
case. This kind of study for other geometries and bigger clusters are under way. 
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Figure 4. K-edge z-polarized absorption cross section for the a state of the Sei 
molecule, for various Z-values up to Z max = 60, calculated by full inversion of the MS 
matrix (I-GT). 



5. Applications 

Application of the present FP-MS theory to two cases which, according to our 
experience, need significant non-MT corrections for a good reproduction of the 
absorption data (i.e. diatomic-linear molecules and tetrahedrally coordinated 
compounds) have already been presented in Ref. [22] for the i^-edge of Se 2 and the Si 
L 2 ,3 edge of crystalline Si02 (a-quartz). There it was shown that a good description of 
the anisotropies of the potential leads to a substantial improvement of the calculated 
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absorption signal in comparison with the experimental spectra. 

In this section we present another application to the A'-edge absorption of Br 2 and 
discuss a preliminary application of the NMT approach to the study of the performance 
of two effective optical potentials, the Hedin-Lundqvist (HL) potential and the Dirac- 
Hara (DH) in the case of a transition metal. 

It should be emphasized that all potentials used here and in Ref. |22j are non-self- 
consistent, since the starting charge density is obtained by mere superposition of atomic 
densities. Therefore the agreement or disagreement with experiments might change if a 
self-consistent charge density were used, although from our experience the effect of this 
latter has a minor impact on the spectra than the elimination of the MT approximation. 
In any case, one of the motivations for pursuing the FP-MS method was exactly the 
study of the performance of the various models of optical potential together with the 
effect of the self-consistent charge density, once that the geometrical approximation of 
the potential had been eliminated. The application of the present real space theory to 
the generation of the self-consistent ground state density using the well known technique 
of contour integration in the complex energy plane is under way. 

In order to obtain the absorption spectra we start from the well known expression 
of the absorption cross section in terms of the GF, given by 

o~tot{w) = —8iTahu! 

x E ^ f (C( r )l^-r|G(r,r / ;£;)|£.r / |«/»L(r / ))drdr' (122) 

For more details and other spectroscopies we refer the reader to Ref. [35]. We used 
all three forms of GF given by Eqs. (93), (94) and (95). While the last two are 
numerically stable and give almost coincident spectra, the first one shows occasionally 
small but noticeable kinks in the calculated spectrum and sometimes small deviations 
around maxima and/or minima of the cross section compared to the other two. This 
is a known phenomenon which is now exalted compared to the MT case, where it was 
almost unnoticeable. It is due the fact that the singularities of the S-matrix in the 
definition of the scattering basis functions $(r) in Eq. (70) and those of T _1 in the 
inverted MS matrix r = (T -1 — G) do not compensate exactly. Therefore, even though 
the three forms are formally equivalent, from a computational point of view, form Eq. 
(93) is to be avoided. 

Fig. (j5|, shows the experimental unpolarized K-edge absorption cross section of 
the diatomic molecule Br 2 [2] in comparison with a NMT and a MT calculation as a 
function of the photo-electron kinetic energy E referred to Eq, the true zero of the non- 
self-consistent molecular potential at infinity. All spectra were normalized at a common 
energy point between 20 and 30 eV. For the NMT case we partitioned the space with 
24 Voronoi polyhedra arranged on a BCC lattice: two of them around the physical 
atoms and 22 empty cells (EC) to cover the rest of the space where the density (and the 
potential) are significantly different from zero. Z max was taken equal to 4 in all polyhedra. 
We gave a small finite imaginary part to the energy of the order of (~ 0.02 eV) in order 
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to be able to use the same Green's function expression for the cross section Eq. (122) 
both for bound and continuum states. To calculate the absorption spectrum, we used 
the real part of an Hedin-Lundqvist (HL) potential and then convoluted the result with 
a Lorentzian whose width is equal to the that of the core hole (2.52 eV). We see that 
the agreement with experiment is rather good. In contrast, the MT approximation of 
the potential turns out to be rather poor. 



Br 2 Gas(300) Exp. 
NON-MT Calc. 
MT Calc. 









10 







10 20 
E - E (eV) 
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40 



Figure 5. K-edge unpolarized absorption cross section for Br 2 molecule, showing the 
comparison between the MT and FP-MS calculations against the experimental data 



We then present in Fig. ([6]) a preliminary application of the NMT approach to 
assess the performance of the HL against a DH potential, assuming that the losses are 
sufficiently well described in both cases by the imaginary part of the HL self-energy, 
in the case of HCP Co metal. As is well known, the real part of the HL potential is 
composed of two terms: the static Hartree-Fock (HF) exchange, known also as Dirac- 
Hara (DH) exchange, coming from the constant part of the dielectric function and the 
dynamically screened exchange-correlation contribution (HLXC), originating from the 
a»-dependent part (see Appendix A of Ref. [4"2]). 

This calculation (and other similar along the same line) were performed without 
any adjustable parameter. In all cases the number of atoms forming the cluster is about 
140-150, lying inside a sphere of about 7-8 A, enough to obtain spectral convergence 
in the presence of the complex part of the potential. The charge density was obtained 
by superposition of neutral atom charge densities, from which the Coulomb and the 
exchange-correlation potential are calculated. In the case of close-packed structure, this 
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E-E (eV) E-E (eV) 

Figure 6. Comparison between Co K-edge absorption calculated with complex HL 
(left) and DH (right) potentials with experimental results. (Color online) 

fact should not be an handicap. 

By contour integral of the Green's Function over the energy range of the valence 
sates, the Fermi Energy was determined to be around -10 eV with respect to E = 0, i.e. 
the zero of the cluster potential at infinity. It serves to define the local momentum of the 
photo-electron in the calculation of the HL (DH) potential, but the calculated spectra 
are rather insensitive to small variations of this quantity by 1-2 eV. No self-consistence 
loop was attempted to find a self-consistent charge. The core hole width was taken into 
account by adding 0.7 eV to the complex part of the potential. 

Surprisingly enough, the comparison shows that the DH potential gives overall 
better agreement with the experiments than the HL one. A similar situation is found 
for other transition metals and has been reported elsewhere. [13] Notice that the same 
conclusion was drawn in Ref. [H] for Cu2MnM, where M = Al, Sn, In, although in 
the MT approximation. 

6. Conclusions 

We have developed a FP-MS scheme which is a straightforward generalization of the 
usual theory with MT potentials and implemented the code to calculate cross sections 
for several spectroscopies, like absorption, photo-electron diffraction and anomalous 
scattering, as well as bound states, by a simple analytical continuation. The key point 
in this approach is the generation of the cell solutions ( r ) for a general truncated 
potential free of the well known convergence problems of AM expansion together with 
an alternative derivation of the MSE which allows us to treat the matrices S and E as 
square, with only one truncation parameter, given by the classical relation Z max ~ k R b . 
The fact that the theory can work with square S and E matrices is of the utmost 
importance, since this feature allows the definition of the cell T matrix and its inverse, 
recuperating in such way the possibility to define the Green's function and to treat a host 
of problems, ranging from solids with reduced symmetry to randomly disordered alloys in 
the context of the CPA, as mentioned in the introduction. In this way one can also show 
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that the wave function and the Green's function approach provide the same expression 
for the absorption cross section for continuum states and real potentials, through the 
application of the generalized optical theorem (see Appendix E). For transitions to 
bound states the two methods are not equivalent, due to the different normalization of 
continuum and bound states, unless one normalizes to one the wave function for these 
latter. However this procedure, although feasible, is rather cumbersome (this was one 
of the reasons for abandoning the MS method in favor of the simpler linearized methods 
in band structure calculations). Instead, the Green's function expression for the cross 
section Eq. (122) can always be used, since it gives the correct normalization in both 
cases simply by analytical continuation. We have exploited this fact when calculating 
the cross section for the Se^ and Br 2 diatomic molecules. 

Moreover, in the present paper we have been able to show that the FP-MST 
converges absolutely in the Z max — > 00 limit (modulo a slight modification of the free 
propagator matrix G which is practically unnecessary) in the sense that the scattering 
path operator of the theory can be found in terms of an absolutely convergent procedure 
in this limit. We have thus given a firm ground to its use as a viable method for electronic 
structure calculation and at the same time have provided a straightforward extension of 
MST in the Muffin-Tin (MT) approximation for the calculation of x-ray spectroscopies. 
Also Quantum Chemistry calculations might benefit from this method in that it avoids 
the use of basis functions sets. 

Finally it is worth mentioning that in giving a new scheme to generate local 
basis functions for truncated potential cells, we have provided an efficient and fast 
method for solving numerically a partial differential equation of the elliptic type in 
polar coordinates, which can also be used to solve the Poisson equation in the whole 
space by the partitioning method. 
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Appendix A. The Mathieu functions 

For the convenience of the reader we give here a brief account of the Mathieu functions. 
The solution of the 3-dimensional Mathieu's equation (24) of the text is obtained by 
separation of variables 

V, z ) = fx(x)f y {y)f z (z). (A.l) 
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in terms of functions / solutions of the one-dimensional Mathieu's equation [27] 
cPf(r) 



dr 2 



(—a + 2q cos 2r)/(r) 



A solution of Eq. ( |A.2| ) having period ir or 2ir is of the form, 

oo 

f(r) = (A m cos mr + B m sin mr) 



(A.2) 



(A.3) 



m=0 



where B can be taken as zero. If the above expression is substituted into Eq. ( A.2 ) 
one obtains 

oo 

2j [(a - m 2 )A m - q(A m _ 2 + A m+2 )] cos mr 

m=-2 
00 

+ 22 t( a ~~ m2 ) B m - q(B m -2 + 5 m +2)] sin mr = 



(A.4) 



m=— 1 



with A_ m 
types, 



if m > 0. Eq. ( A.4 ) can be reduced to one of four simpler 

(A.5) 



fo( r ) = A ?™+p cos ( 2m + P) r i V = or 1 

m=0 
00 

fi(r) = 2J ^2m+ P sin(2m + p)r, p = or 1. 



(A.6) 



m=0 



If p = 0, the solution is of period 7r; if p = 1, the solution is of period 2n. f Q is an even 
solution, and fi is an odd solution. The recurrence relations among the coefficients of 
these basic solutions are easily obtained from the general relations Eq. ( A.4 ). For even 
solutions of period n we find 



aA — qA 2 = 

(a - A)A 2 - q(2A + A 4 ) = 

(a - m 2 )A m - q{A m _ 2 + A m+2 ) = 0, m > 3 

and of period 27r, 

( a _l)^ 1 _ g (A 1 + A 3 ) = 

(a - m 2 )A m - q(A m _ 2 + A m+2 ) = 0, m > 3. 
For odd solutions of period ir, 

(a - 4)B 2 - gA 4 = 

(a - m 2 )B m - q(B m _ 2 + £ m+2 ) =0, m > 3 

whereas for period 27r, 

(a - 1)5! + q(B 1 - B 3 ) = 

(a - m 2 )B m - q(B m - 2 + =0, m > 3 



(A.7) 
(A.8) 
(A.9) 

(A.10) 
(All) 

(A.12) 
(A.13) 

(A. 14) 
(A.15) 
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Table Al. First few Eigenvalues of Mathieu functions for different q values 



parity 


even 


odd 


period 


7T 


2tt 


7T 


2tt 


q=0.01 


a =-4.99995xlCT 6 
a 2 = 4.00004 


ai = 1.00999 


6 2 = 3.99999 


6i = 0.989988 


q=0.02 


a =-l-99991xl0- 5 
a 2 = 4.00017 


oi = 1.01995 


b 2 = 3.99997 


by = 0.97995 


q=0.03 


a 0= _4.49956xl0- 5 
a 2 = 4.00037 


oi = 1.02989 


6 2 = 3.99993 


by = 0.969888 


q=0.04 


a =-7.9986xl0- 4 
a 2 = 4.00067 


ay = 1.0398 


b 2 = 3.99987 


by = 0.959801 


q=0.05 


a =-1.24966xl0" 3 
a 2 = 4.00104 


oi = 1.04969 


6 2 = 3.99979 


by = 0.949689 


q=0.1 


a =-4.99454xl0" 3 
a 2 = 4.00416 


oi = 1.09873 


b 2 = 3.99917 


by = 0.898766 


q=0.2 


a =-l-99133xl0- 2 
a 2 = 4.01658 


oi = 1.19487 


b 2 = 3.99667 


by = 0.795124 


q=0.3 


a =-4.4566xl0" 2 
a 2 = 4.03706 


ai = 1.28832 


6 2 = 3.9925 


by = 0.689166 


q=l 


a =-0.455139 
a 2 = 4.3713 


oi = 1.85911 


b 2 = 3.91702 


by = -0.110249 


q=2 


a =-1.51396 
a 2 = 5.17267 


oi = 2.3792 


6 2 = 3.67223 


by = -1.39068 


q=5 


a =-5.80005 
a 2 = 7.44911 


oi = 1.85819 


b 2 = 2.09946 


by = -5.79008 


q=10 


a =-13.937 
a 2 = 7.71737 


oi = -2.39914 


b 2 = -2.38216 


by = -13.9366 



It is convenient to separate the characteristic values a into two major subsets: 

a = a r , associated with even periodic solutions 
a = b r , associated with odd periodic solutions 



where r describes the index of the eigenstate. Table XT gives the first three eigenvalues 
associated to even periodic solutions and the first two associated to odd periodic 
solutions (bo = 0), for some selected values of q. They can serve to generate the 
Mathieu functions using the above recurrence relations to determine the coefficients 



in the solutions (A.5) and (A. 6) 



Full- Potential Multiple Scattering Theory with Space-Filling Cells for bound and continuum states39 
Appendix B. Asymptotic behavior of KKR Structure Factors 

For v — > oo through real positive numbers (in practice for v ^> |z|), the other variables 
being fixed, one has [27] 



2nuJ \2u 

-iHt-^i^T (b.i; 



XttuJ \2v 
where e is the Neper number. Remembering that 

jn{z) = J— J n+1/2 (z)\ 



h±(z) = ] [^H± +l/2 (z) (B.2) 



we find for the asymptotic behavior of the spherical Bessel and Hankel functions 

z n , / 1 \ "+ 1 



We need to find un upper limit for G l [ L , given by 

G% = —Airik ^ V+l " C{L, L'; L") h+(p) Y v ,(p) (B.4) 

L" 

where p = kRij, p = and C(L,L';L") are the Gaunt coefficients. To establish an 
upper limit for this expression when I is fixed and I' > p we replace each hfi, (p) | in the 
sum by its maximum value |^tu/(p)|, use the asymptotic value in Eq. (B.3) and the 
relation £ L „ C(L, V; L") Y L „(p) = Y L (p)Y L ,(p) to obtain 



\G%\ <4nk\ht +l ,(p)\ J2\C(L,L';L")Y L „(p)\ 

L" 

« Ank I ht +v (p) \\J2C(L,L'; L") Y L „ (p) | 

L" 

= 4nk\ht +l ,(p)\\Y L (p)Y L/ (p)\ 

< [(2/ + l){2l> + l)] 1/2 p ^ +1 e J +1/2 W + 1') + (B.5) 

since \Y L (p)\ < y/(2l + l)/(4vr). Notice that the approximation ^2 L „ \C(L, L'; L")Y L n(p)\ 
\^2 L „ C(L,L, L")Yl»(p)\ entails only errors 0(1) in all I variables, as can be ver- 
ified by explicit calculation, and therefore completely negligible with respect to the 
power behavior of the rest of the factors. In any case, since \C(L, L'\ L")Y L n(p)\ < 
[(2/ + 1)(2/ , + 1)] 1 / 2 /(4tt) El»(2Z" + 1), at the cost of introducing a non influential extra 
factor [2(1 + I') + l] 2 in Eq. (B.5) we would get a rigorous inequality. This expression is 
obviously also valid for / ^> p. 
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Under the same conditions, assuming I' ^> I we derive 
\J l l L> \<^\jv-i{p)\\Y L {p)Yv{p)\ 
< [(2/ + l)(2/' + l)] ] 



j'-i 

11/2 P I' -1+1/2 



V2 



[2(i' -i) + iy- i+i 



(B.6) 



The inequalities Eqs. (B.5 ) and (B.6 ) can be used to obtain other useful inequalities 
used throughout the paper. For example, for fixed I, using again Eq. (B.3), one obtains 

k 



\G&,M*i)\ < 



Rij 



{klk. 



(B.7) 



implying that the series H^(vi) = ^2 L , G % l L ,Jy{jj) is absolutely and uniformly 
convergent in the angular domain. The uniform convergence comes from the application 
of the Weierstrass criterion (see Ref. [ID], sect. 3.34, pag. 49). 
Similarly one finds 

I'+i 



\J? L ,H+(v )\ < 



Rir 



k 



(kR t 



u+i W - + 1]' / 



21 + 1 



47T 



(B.8) 



showing that the series H^iv^ k) = J2l> ^ll'^l'^o] k ) is a l so absolutely and uniformly 
convergent if (r > Rio)- 

Along the same lines we can estimate an upper bound for the atomic T-matrix for 
/, V 3> kR b . We find from Eq. (42) to first order 

"Rb 



J v (v)V(v)^j l (v) d 6 i 



r 2 j v {kr)V vll {r)R LIIL {r)di 

fR b 

J2C(L',L;L") / r 2 jV(Ar)M0.H(*r)<k 

TP JO 



(B.9) 



where the last step follows from the fact that under the above assumptions Rl'l ~ ji&LL 1 - 
Taking into account that C(L',L;L") ~ l/\/47rO(l) for all L-values and using again 
Eq. (B.3) we obtain 

\T LL A <4ZZ' 



Rb 



All'k 



r 2 \ji'(kr) \ \ v \i-i'\( r ) \ \ji(kr)\ dr 
j+i'+i 



(21 + (21' + 



Rb 



x / r l+v+2 \V^ V \(r)\ dr 



< 



J eff 



Ml' 



(kRi 



\l+l'+2 J+V+l 



k 2 i + v + 2 (2i + iy+ i (2i' + 1Y+ 1 ^ B ' 10 ^ 

with the understanding that Vi = Viq, assuming that |V/(r)| < 2Z e ff/r in atomic units 
and that |V/(r)| is decreasing with I. 
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Based on the above inequalities we easily obtain 

(ll'fl 2 



G%T LV \<%V2e l ' 2 Z eff R, 



ffIXb l + V + 2 
R b \ l+V+1 (21 + 21' + 

Wj) (2/ + (2V + 

Specializing to the case where / is fixed and I' is running, we also find 



r (ll') 1/2 l' 
\G 1 I l ,T l >l>\ < 4 Z eff 



x 



Z' + l 



D 2l'+l 



(21 + 21' + I 



d+l' 



(B.ll) 



(B.12) 



(kRijY +l ' +1 e l + l '+ 1 / 2 (2V + l) 2 ^ 1 
which is useful in discussing questions related to the convergence of MST. 

Finally we note that all the above inequalities and convergence conditions remain 
valid for complex arguments p, provided it is replaced by its module \p\. 



Appendix C. Surface identity for scattering states 

In the case of short range potentials (i.e. potentials that behave like l/r 1+<E with positive 
e as r — > oo) the Lippmann-Schwinger equation for scattering states at energy E = k 2 

</>(r; k) = o (r; k) + J dr'G (r - r'; k)V(v')^(v' ; k) (C.l) 

is a consequence of the Schrodinger equation 

(V 2 + E- V(r))ip(r;k) = (C.2) 
together with the relations (0o( r ; k) = e 



(V 2 + £)<Mr;k) = 

(V 2 + E)G (r - r'; k) = 5(r - r') 



Starting from Eq. (C.l), we derive the identity 



/ dv' [G {r - r'; k)V(v') - 5(v - v')\ ^(r' ; k) = -0 o (r; k) 
Jn 



(C.3) 
(C.4) 

(C.5) 



where f2 indicates the whole space. Using Eq. (C.4) to replace the delta function, and 



the Schrodinger equation (C.2) to eliminate V^(r') we obtain 

N+l „ 

V / {G (r-r';fc)(V 2 + i?Mr';k) 
j= i J Vj 

- ^(r'; k)(V 2 + E)G (r - r'; k)^ = -0 o (r; k) 
where we have decomposed the whole space as Q = ^2^=1 sucn that Qn+i = 

Transforming to surface integrals by application of the Green's theorem 

N+l 



f\ [ [G (r-r';k)Vij(r';k) 

3=1 Js i 

-^(r';k)VG (r - r' ; k)} ■ n'^ 



(r;k) 



(C.6) 
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We now observe that the surface integral over the surface Sjv+i of the volume 
Q N+1 = Q Q has two contributions, one coming from the surface S Q of X^Li the other 
one at infinity, as the limit as R — )■ oo over the surface of a sphere S^, of radius i?. 
This latter is easily calculated on the basis of the asymptotic behavior of ip(r; k) in Eq. 
(41) and the expansion (36) and gives exactly — </>o(r; k), canceling the rhs term in Eq. 



(C.6). Therefore we recover the identity (32) of Section 3.1 

N 

V / {C (r-r',A ; )V^(r / ;k) 

— ip(r'; k)VGo(r — r'; k)} ■ n'jda'j 
= [ {G (r-r',fc)VV;(r';k) 

J So 

- ip(r'] k)VG (r - r'; k)} ■ n' da' 



Appendix D. The Generalized Optical Theorem 



(C.7) 



For convenience of the reader we give here a proof of Eq. ( 79 ) in the case where T = 



i.e. when empty cells cover the volume f2 Q up to the point at which the asymptotic 

(D.l) 



behavior in Eq. (41) begins to be valid. We start by observing that 

k 



7T 



dkP L ,(k) [I£(k) 
so that, using the relation Eq. ( J77[ ), we find 

dksi,(k) [si(k)]* = 5>^,j^ 



'71% 

r A'L' 



LA 



k 
7i 



(D.2) 



where we have used the symmetry of r. Based on the relations Eqs. (102), (103) and 



(]105|), valid at any energy, and due to the reality of the matrices K, N and J for real 

(D.3) 



potential, we can write 

T = kT l [K -N + iJ]' 1 



so that the rhs of Eq. (D.2) becomes 
k , _,,■„■ 11 



7T 



{tJtY 3 



LL' 



7T 2l 



LL' 



7T 



LL' 



(D.4) 



in keeping with Eq. ( 79 ) . 



Appendix E. Wave function and GF equivalence for absorption cross section 

In the independent electron approximation, the core level photoelectron diffraction 
(PED) cross-section for the ejection of a photoelectron along the direction k and energy 
E = k 2 from an atom situated at site i is given by [35] 
da 



dk 



8Tr 2 ahuJ2\(®^( r i-M£-ri\<i>lS r i)}\ 



(E.1) 
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Here is the time-reversal operator, e the polarization of the incident photon and 
<j) c Lc (ri) the initial core state of angular momentum L c (we neglect for simplicity the 
spin-orbit coupling, which can be easily taken into account). Due to the localization of 
the core state, we need only the expression of the continuum scattering state in the cell 
of the photoabsorber, given by 



^(r l ;k) = ^ J Bl(k)$ L (r i 



so that 



da 
dk 



\-n 2 a 



(E.2) 



(E.3) 



where B l L (k) is given by Eq. (77) and we have defined the atomic transition matrix 
element 



M Lc l{E) 



dr<p c Lc (r)i -r$ L (r) 



(E.4) 



The total absorption cross-section, in the case of real potentials, is obtained by 
integrating the PED cross-section over all directions of photoemission 

2 



J dk^Z = 8n 2 ah-uY, J dk 



(E.5) 



m c LL' 



by application of the optical theorem (79). This is exactly the form that one would 
obtain starting from Eq. (122) and using the expression (93) for the GF. 



Appendix F. Exploitation of point symmetry 

In a cluster (to which we shall also refer as a molecule), point symmetry can be 
used to advantage to simplify the problem and reduce the size of the MS matrix. 
Specifically we consider the case where the cluster remains invariant under a finite 
group of transformations Q relative to the molecular center R Q . This group has a finite 
number of finite-dimensional irreducible representations (irreps) Tj (j = 1, 2, g). Due 
to the symmetry, the cluster will consist of V groups of equivalent atoms, transforming 
into one another under the operations of the group; and for each group p — 1,2, V, 
there are N p atoms labeled by i p . If iV is the total number of atoms in the cluster, then 

Under these assumptions there exists a unitary transformation C that block- 
diagonalizes the MS matrix according to the irreps. For each value of the angular 
momentum index I, this transformation is labeled by the angular projection m and the 
site i p on one side, and on the other the irrep, the row p of the irrep and a further 
index n that distinguishes independent orthogonal symmetrized basis functions with 
the same /. The matrix elements of C are easily obtained by applying the projection 
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p 

operator J2rM P p(R) Pr to a spherical harmonic function Yi m (r ip ) 9(Rb — r ip ) centered 
on the site i p and defined on the surface of the bounding sphere Rb of the cell fL Here 
Pr is the generic operation belonging to the group Q corresponding to the coordinate 

p . 

transformation R, and M P p (R) is the matrix element corresponding to R in the matrix 
representation of irrep Tj of the group, p labeling the row of the irrep. As usual in group 
theory [15] the effect of Pr on the function /(r) is given by /(-R -1 r). In this way, if the 
result is not zero, one generates a symmetrized spherical harmonic function given by 

Ki!' P (r p )= Y,M t p ^R)PrYUK) 

R 

= £ M%(R)J2 D U(R)YUii' p ) 

R fi 

m,i p 

where D l (R) is the Wigner rotation matrix corresponding to the transformation R jl5] 
and i' p = Ri p . Due to the orthogonality of the basis functions 

/«„(?,) H„.(f,,)dn =s mm ,s v , 

r p' 

J <i' P ( r p )K V Z;\ r p ) dn = 5 IV 5 nn , 8 Tj r., <V ( F - 2 ) 
we obtain CC = CC = I, i.e. 

E tfn,'m C l'n',m = S W <W b~TV <5 V 

m, i p 

E E C in,m C ln,'m> = S rnm> Si p i> p (F.3) 
Tp n 

where for simplicity we have dropped the index j from the symbol T of the irreps. 

Now, if M = M % l L i = {T~ 1 —G) l [ L i is the MS matrix in the non-symmetrized site and 
angular momentum indices, its symmetrized version is given by M s = C MC. Therefore 
for any representation T we have, putting for short A = (I, n) and remembering that 
L = (Z,m), 

ip m, m' 

Np Nq 

^r,pg _ \ ^ \ ^ \ ^ ^r,i p n ipi q n T,i q /-p r\ 

"A, A' Z^Z^ °A,£ L, V A', V \ t - b ) 

ip iq m, m! 

Here Tj ,p describes total scattering power of group p and G v s ' pq is the symmetrized 
matrix of the KKR structure factors. The presence of an outer sphere contribution JT° J 
is treated on the same footing and contributes C JT J C = C JCCT CC J C = J S T S J S 
for each representation V. Since the outer sphere is centered at the origin of the cluster, 
it has no partner spheres equivalent to itself. 
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All these matrices are labeled only by the groups of equivalent atoms (prototypical 
atoms), the angular momentum I and possibly the index n mentioned above, realizing a 
sizable reduction in dimensions. Notice that, since the molecular hamiltonian is invariant 
under the operations of Q, the symmetrized matrix elements do not depend on the row 
p of the representation V. Moreover, in order to find the symmetrized T-matrix relative 
to an equivalent group of atoms, we do not need to calculate T^ L , for all sites in the 
group, since these are related to one another by the relation 

T% = E D U R ) T ii>»'D l m ,^(R) (F.6) 
fin' 

where, as before, i' p = Ri p . This relation is a consequence of the invariance of the 
potential and the T-matrix T tp (i,r') under the operations of the group Q. Specifically 

V£ v = f r L (?)V^(r,r)y L ,(r)dr 
Y L {r) V^(r,.Rr)y L /(r)dr 
= J V^iT 1 ?) V^^y^iT^dr 

= Dlm(R)V l % fl ,D i ;, m ,(R) 

valid also for the matrix elements T l £ L , = J Yl(t) T tp (i, r') Yli(j') d r dr'. Notice that 
the transformation and symmetrization properties are the same for T and T" 1 , so we 
can act directly on this latter. 

In the MT case the T-matrices are angular momentum diagonal and m and site 
independent within a set of equivalent atoms, so that 

rpV,p _ rpp \ ^ \ ^ ^r,i p ^F^p 

1 A,A'— 1 l / , / ; °A,L °A',L 

i p m 

= Tf5 KA , (F.7) 
The group T p -matrix can also be calculated directly from the symmetrization of the 



radial part of the basis function Rll'^^), which transforms as T in Eqs. (F.4) and 



(F.6). Even though the number and type of operations to perform are exactly the 
same as for obtaining T p , in this case there is the added advantage of generating a 
symmetrized form of the scattering wave function, needed for example to calculate the 
absorption or photo-emission cross-section. 

The symmetrization of the local wave function in a group of equivalent atoms for 
an irrep T is obtained by observing that the following function ^(r^) is invariant under 
any operation of the group 

^p) = E E A L R L'L(n p ) Y L ,(T ip ) 

= (A\R\Y) 

= (A\CCRCC\Y) 
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= Y^ A A p &*U r p) K % p (*p) (f.8) 

whereby 



AA' 



An 



-^AA'( r p) — E E ^A,L Rl,L> ( r p) ^A',L' 



i p m, m' 



Inside the MT sphere -X" A , A (r p ) = r -^A'A( r p) * s solution of the symmetrized equation 



(14) 



E 

A" 



~ Z(Z + 1) . . 



X^ A ,(r) 



(F.9) 



where we have written for simplicity r for r p and V^ A „(r) is given by 



An 



(F.10) 



«d m, m' 



Ti being the identical representation of the group Q. Near the origin X^£,, (r) ~ 
rji(kr)Kl' p (r p )5 AA >. 

Across the truncated boundary, we instead use the symmetrized version ofEq. Q, 
so that putting P A ' p (r p ) = r$ A p (r p ) and dropping again the index p 

^ " 1 " ^M(r,*) (F.ll) 



<ir 2 



+ £-F(r,r) 



P A r (r,r) 



where 



L 2 P A r (r, f) = Y, HI' + l)rR r A , A (r)KUi [ 



(F.12) 



and we use starting values given by Eq. (F.9). Equation (F.ll) is obtained from Eq. 
([!]) by applying on the left the projection operator ^2 R Mpp (R) P R , taking into account 
thatV(r tp ) = V(R- 1 r l ,). 

In terms of i? AA , (r) is then possible to define the symmetrized version of the matrices 
E % l L , and S l [ L , as 

E P AAI = (R P b ) 2 W[-i K h+,Ri AI ] (F.13) 
Si A > = (R P b ) 2 W[j h Rl A ,} (F.14) 
and derive the symmetrized equivalent of all the quantities introduced toward the end 



of section (3.1). In particular the amplitudes P A (k) are solutions of the symmetrized 
MSE 



E l^A'A' + G l P A B l^) = J A( k ) 



(F.15) 



qA> 
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where I^(k) = ^2 mip C^'m -^mOf)- Assuming that the photo-absorber is located in cell 
Q Q at the origin of the coordinates R D , the symmetrized PED cross section for the final 
state irrep T with degeneracy d(T) takes the form 



da v 
dk 



where the atomic dipole transition matrix element 



Mlf A (E) 



fc(r) L> r «(r) <^(r)dr 



(F.16) 



(F.17) 



obeys the selection rules of the Wigner-Eckart theorem jlS] for the finite group Q. We 
have assumed that the dipole operator transforms according to the irrep Yd- 

Appendix G. Finiteness of Tr (K\K 8 ) 



In this appendix we show that Tr {K\K^) is finite. Starting from Eq. (114), we partition 



the space and the potential in the way described at the beginning of Section |3.1| and 
define a new kernel K that coincides with the kernel K s for r and r' in different cells 
and vanishes identically when r and r' happen to be in the same cell. 

For this new kernel, Tr(K^K) is finite (< N < oo) for an even larger class of 



potentials than that defined by Eq. (114) so that, rewriting K(r, r') in operator notation, 
we have 



N > j dr (rlK^Klr^ 

= I I I dEdE'dE" [ dr ( v \ J l(E)) 

J J J LL'L" ^ 

x (j L {E)\tf\J L ,{E')) (j L ,(E')\K\J L „{E")) (J L »{E")\x) 
dEdE'J2\^Lu(E,E')\ 2 

LL' 

> J2\Kll>(E,E)\ 2 (G.l) 




LL' 



taking into account that the functions Jl(E)(t) = (A;/7r) 1/ ' 2 ji(A;r)lL(r), with the 
normalization to one state per Rydberg, form a complete orthonormal set. 

Now it is easy to see that the matrix Kll' is the asymptotic form of the kernel K s 



in Eq. (118) for high values of the indices LL', since 
K W {E,E) 

= V / dr t J L (E)( ri )\V t (^ ll/2 * 



VAYi 



I diJG+(r, - r'j + R«; kW^y^v'^E)^) 

f ^ME^rm^M^um^Gfv 



ijtj AA' 
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x 



-EEKa] 1/2 ^aa' [ T ^] 1 ' 2 (G.2) 

i^j AA' 

the last line following by the fact that asymptotically, when LA,L'A' 3> kR\ (Vi), 

[T2 A ] 1/2 ~ f d^ME^rmi^^UE)^) (G.3) 

Notice also that we have used the two center expansion for the free Green's Function 

G+(r - r'; k) = J A (£)(r 4 ) G l [ AI J A (£) (G.4) 

which might diverge if 7*3 + rj > Rij, e. g. for neighboring cells. A similar problem 
is encountered when formulating the variational derivation of MST (see, for example, 
section 6.5.3, page 140 and following of Ref. [13]). One way to solve it is to use the 
displaced cell approach [13J, whereby one can write 

Kll>(E,E) 

= EE|EKa] 1/2 JAA(b)G AA ,(R^+b) [n, L ,] 1/2 | (G.5) 

i¥=j A I AA' J 

provided that |Ry + b| > R\ + R J b and the sums inside the curly brakets be performed 



first. Here J AA (b) is the usual translation operator in MST, given by Eq. (51) with the 



vector Rj, replaced by the vector b. The tilde over the symbol in Eq. (G.4) was 



meant to be a reminder to use this procedure. Notice that the vector b depends only 
on the geometry of the partition of the space in cells and is independent on I. 



In this way the expression Eq. (G.5) is always convergent and is such that 



\K LL ,(E,E)\ 2 = Tr(K^K) is finite. Consequently also Tr (KtK s ) is finite. 
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